These  fioures  show  one  of  the  several  remarkably  stable  and  robust  structures  of  the  equations  that  model  double  diffusion 
convection.  In  such  flows,  thermal  effects  and  density  stratification  (salinity)  compete  to  form  complex  equilibriurn  structures. 
The  snapshots  in  this  figure  show  the  cross  section  of  the  state  of  a  periodic  solution  known  as  *dscillating 
the  structure  is  neither  a  standing  wave  nor  a  traveling  one.  For  more  details  of  the  possible  patterns  in  the  double  dinusion 
equations,  see  the  article  by  Y.  Renardy  in  this  issue. 
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Dynamical  Systems  and 
Ocean  Dynamics 


Reza  Malek-Madani,  Guest  Editor 
Office  of  Naval  Research 


Introduction 

Since  the  seminal  works  of  Liapunov  and  Poincare,  re¬ 
search  in  dynamical  systems  has  provided  scientists  with  a  set 
of  practical  tools  to  investigate  the  qualitative  behavior  of 
solutions  of  differential  equations.  In  the  past  fifty  years  this 
area  of  mathematics  has  produced  some  of  the  most  insightful 
and  precise  information  about  the  nature  of  the  dynamics 
present  in  many  important  applications  ranging  from  biology 
to  physics.  Recently,  a  group  of  applied  mathematicians 
joined  together  in  a  concerted  effort  to  develop  dynamical 
systems  tools  appropriate  for  problems  in  ocean  dynamics. 
The  articles  that  follow  will  give  an  introduction  to  these 
methods  and  provide  a  sample  of  the  type  of  questions  that  one 
can  properly  pose  and  answer  in  a  dynamical  system  setting, 
and  point  to  the  future  directions  of  research. 

The  mathematical  modeling  of  the  ocean  is  a  relatively 
young  science  that  still  takes  many  of  its  cues  from  observation 
and  experimentation.  Although  the  basic  governing  equations 
have  been  well-known  for  many  decades,  a  large  class  of 
partial  differential  equations  are  present  in  the  literature  whose 
derivations  depend  on  the  special  scales  in  the  geometry  of  the 
region  of  interest  as  well  as  the  strength  of  the  forces  involved. 
In  addition  to  the  information  one  draws  from  these  models,  a 
substantial  amount  of  data  is  collected,  either  remotely  or 
in-situ,  that  has  contributed  significantly  to  our  understanding 
of  the  basic  phenomena  as  well  as  challenged  some  of  the 
underlying  assumptions  that  have  led  to  the  mathematical 
models.  One  approach  to  reducing  the  gap  between  the  sophis¬ 


tication  and  physical  realism  of  the  mathematical  models  and 
the  information  gathered  in  real  data  is  to  develop  mathemati¬ 
cal  tools  that  are  not  so  heavily  model  dependent  that  they  lose 
their  significance  when  one  mathematical  model  is  discarded 
in  favor  of  another.  With  this  point  in  mind,  the  methods  and 
techniques  outlined  in  the  next  three  articles  in  this  journal  aim 
at  understanding  the  basic  issues  that  describe  transport  and 
mixing  in  the  ocean  and  atmosphere  as  well  as  quantifying  the 
competition  among  the  forces  that  lead  to  double  diffusion 
convection. 

Computational  models  have  been  traditionally  the  widely 
used  tool  for  predicting  the  response  of  the  ocean  and  the 
atmosphere  to  various  mechanical  and  thermodynamical  proc¬ 
esses.  In  many  respects,  these  efforts  have  closely  followed 
the  development  of  the  numerical  analysis  of  the  Navier- 
Stokes  equation  and  have  attempted  to  find  the  appropriate 
range  of  various  viscosity  parameters  that  would  lead  to  rea¬ 
sonable  numerical  simulations.  This  strategy  has  paid  off 
well  in  the  context  of  low  Reynolds  number  dynamics  mod¬ 
eled  by  this  equation  while  it  is  a  rich  area  of  research  for  high 
Reynolds  number  turbulence.  Because  the  Navier-Stokes 
equations  are  generally  considered  the  correct  mathematical 
model  for  many  applications,  it  is  felt  that  the  investment  of 
resources  in  developing  more  accurate  numerical  schemes  for 
this  specific  equation  is  worthwhile.  (See  reference  1  for  a 
viable  alternative  to  a  direct  discretization  of  the  Navier- 
Stokes  equations,  whether  using  finite  difference  methods, 
finite  elements,  or  the  various  quasi-spectral  methods  that  are 
currently  being  used  in  the  literature).  It  is  not  clear  that  this 
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strategy  is  correct,  or  even  desirable  in  ocean  dynamics,  con¬ 
sidering  the  fact  that  more  than  one  candidate  for  the  “stand¬ 
ard"  model  exists.  It  then  seems  reasonable  that  the  solutions 
of  the  governing  differential  equations  are  understood  at  a 
level  beyond  the  usual  listing  and  tabulating  the  various  values 
of  the  parameters  and  instead  search  for  the  coherent  structures 
that  are  commonly  present  in  these  models. 

This  is  not  to  diminish  the  role  of  numerical  analysis  in 
ocean  dynamics.  In  fact,  quite  the  opposite  result  is  antici¬ 
pated.  The  intent  of  the  new  dynamical  systems  methodology 
is  to  elucidate  the  common  complex  structures  in  the  set  of 
relevant  systems  of  partial  differential  equations.  Once  the 
coarse  properties  of  these  structures  are  identified,  numerical 
simulations  on  a  reasonably  small  and  computationally  tracta¬ 
ble  number  of  initial  conditions  will  reveal  the  finer  properties 
of  the  mathematical  models. 

Scope  of  Applications 

The  mathematics  described  in  the  following  articles  is 
divided  in  two  categories,  depending  on  whether  the  models 
have  a  Lagrangian  formulation  or  an  Eulerian  one.  In  the 
articles  of  S.  Wiggins^  and  C.  Jones  et.  al?  of  this  issue,  the 
goal  is  to  understand  the  nature  of  the  solutions  of  the  per¬ 
turbed  systems  of  ordinary  differential  equations  that  arise  in 
the  Lagrangian  setting  of  motions  in  the  ocean  and  atmos¬ 
phere.  This  point  of  view  is  specially  significant  in  ocean 
dynamics  in  light  of  the  type  of  data  taken  from  the  so-called 
Lagrangian  drifters  that  provide  real  data  about  various  fluid 
particles  in  the  ocean.  Because  only  a  finite  (and  a  relatively 
small)  number  of  these  drifters  are  available  at  the  present 
time,  one  is  left  with  the  challenging  task  of  gathering  infor¬ 
mation  about  a  portion  of  the  ocean  from  a  handful  of  particle 
trajectories.  This  fact  alone  shows  the  acute  need  for  a  kind 
of  mathematical  analysis  that  goes  well  beyond  the  tabulation 
of  parameter  values. 

The  methods  that  are  discussed  in  references  2  and  3  show 
how  certain  structures  associated  with  invariant  regions  of  an 
unperturbed  dynamical  system  are  transported  by  the  per¬ 
turbed  velocity  fields.  Specific  and  delicate  mathematical 
tools  are  developed  to  monitor  the  transport  of  information  for 
a  rather  large  class  of  perturbations.  One  of  the  intriguing 
problems  that  will  be  addressed  in  the  future  by  these  investi¬ 
gators  is  what  happens  to  these  dynamical  systems  methods 
when,  as  is  the  case  in  real  life,  the  velocity  fields  that  the 
differential  equations  are  based  on  are  only  known  on  a  dis¬ 
crete  set  of  points  in  the  ocean.  A  related  problem  is  whether 
one  can  develop  systematic  mathematical  tools  that  can  dis¬ 
tinguish  the  effect  of  chaos  from  stochasticity  in  data  sets  that 
exhibit  complex  and  coherent  structures. 

The  article  by  Y.  Renardy  of  this  issue'^  addresses  the  set 
of  partial  differential  equations  that  model  double  diffusion 
convection.  This  paper  shows  how  to  use  the  center  manifold 


theorem  of  dynamical  systems  and  classify  the  stability  prop¬ 
erties  of  certain  important  solutions  of  the  original  system.  A 
central  feature  of  this  theorem  is  in  its  ability  to  reduce  the 
analysis  of  the  dynamics  to  a  lower  order  system.  Various 
traveling  wave  solutions  of  the  original  governing  system  of 
partial  differential  equations  are  identified  and  numerically 
simulated. 
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Lagirangian  llranspoirt  in 
Geophysical  Flows: 

New  i^proaches  and 
Techniques  from 
dynamical  Systems  Theoiy 

Stephen  Wiggins 

Applied  Mechanics,  and  Control  and  Dynamical  Systems 
California  Institute  of  Technology 


Introduction 

The  standard  map  is  a  two-dimensional  area-preserving 
transformation,  or  dynamical  system,  that  has  been  the  subject 
of  literally  hundreds  of  papers  in  the  dynamical  systems  com¬ 
munity  over  the  past  twenty  years.  The  map  is  given  as  follows 

X  x-fy-TT, 

y  ^y+A:sin(x+y), 

and  it  is  periodic  in  x  and  y.  In  Figure  1  we  plot  a  single  orbit 
of  this  map  for  different  initial  conditions  in  each  of  the  four 
panels  for  the  parameter  value  k=^\.  (By  the  term  “orbit”  we 
mean  a  sequence  of  points  obtained  by  iterating  a  given  initial 
point  according  to  the  map.  Thus  an  orbit  of  a  map  is  a  discrete 
time  trajectory.)  The  orbits  in  the  upper  left  hand  and  lower 
right  hand  panels  explore  a  significant  fraction  of  the  phase 
space.  Moreover,  the  motion  is  chaotic  in  the  sense  that  nearby 
orbits  separate  at  an  exponential  rate.  However,  note  the 
“holes”  in  the  regions.  The  chaotic  regions  do  not  fill  all  of  the 
space.  Note  the  difference  in  the  orbits  in  the  upper  right  hand 
and  lower  left  hand  panels.  These  orbits  appear  to  lie  on  closed 
curves.  There  are  four  discrete  closed  curves  in  the  former  case 
and  one  closed  curve  in  the  latter  case.  Orbits  on  these  curves 


separate  at  a  linear  rate.  These  closed  curves  are  examples  of 
“KAM  tori”  and  are  boundaries  which  orbits  cannot  cross. 
KAM  tori  are  examples  of  lower  dimensional  “phase  space 
structures”  that  evidently  have  a  marked  influence  on  transport 
in  phase  space. 

The  chaotic  regions  are  built  on  the  skeleton  of  a  lower 
dimensional  phase  space  structure.  In  Figure  2  we  show  some 
intertwining  one  dimensional  curves  that  lie  in  the  chaotic 
region  in  the  upper  left  hand  panel  of  Figure  1.  These  curves 
are  the  one  dimensional  stable  and  unstable  manifolds  of  the 
saddle  type  fixed  point  in  that  region.  They  are  the  mechanisnt 
giving  rise  to  the  chaos,  and  they  also  govern  the  transport 
through  such  chaotic  regions  and  play  a  very  important  role  in 
understanding  this  phenomena. 

Thus,  in  this  example  we  see  that  the  phase  space  is  a 
mixture  of  “regular”  and  chaotic  regions,  a  situation  that  we 
now  know  is  typical  for  many  nonlinear  dynamical  systems. 
Moreover,  lower  dimensional  “phase  space  structures”-KAM 
tori,  stable  and  unstable  manifolds  of  periodic  orbils-form  the 
skeleton  in  phase  space  upon  which  the  regions  of  regular  and 
chaotic  motions  exist. 

Despite  the  simplicity  of  this  example,  the  same  mathe¬ 
matical  framework  and  ideas  have  implications  and  uses  for 
many  problems  in  geophysical  fluid  dynamics,  and  we  next 
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want  to  motivate  this  setting  and  approach.  Suppose  one  is 
interested  in  the  motion  of  a  passive  tracer  in  a  fluid  (e.  g.  dye, 
temperature,  or  any  material  that  can  be  considered  as  having 
negligible  effect  on  the  flow),  then,  neglecting  molecular 
diffusion,  the  passive  tracer  follows  fluid  particle  trajectories 
which  are  solutions  of 

i  =  'u(x,r,^),  (1) 

where  t)(jc,r;)i)  is  the  velocity  field  of  the  fluid  flow, 
X  e  IR^,  n  =  2  or  3,and  p  €  IR^  represent  possible  parameters. 
When  viewed  from  the  point  of  view  of  dynamical  systems 
theory,  the  phase  space  of  (1)  is  actually  the  physical  space  in 
which  the  fluid  flow  takes  place.  Evidently,  “structures’’  in  the 
phase  space  of  (1)  should  have  some  influence  on  the  transport 
and  mixing  properties  of  the  fluid.  To  make  this  more  precise, 
let  us  consider  a  fluid  mechanically  more  simplified  situation 
that  can  be  related  directly  to  the  standard  map  example 
described  above.  Suppose  that  the  fluid  is  two-dimensional, 
incompressible,  and  inviscid,  then  we  know  that  the  velocity 
field  can  be  obtained  from  the  derivatives  of  a  scalar  valued 
function,  ,X2,/;p),  known  as  the  stream  function,  as  follows 

Figure  1. 

Four  orbits  of  the  standard  map  for  k=  1. 


In  the  context  of  dynamical  systems  theory,  (2)  is  a  time-de- 
pendent  Hamiltonian  vector  field  where  the  stream  function 
plays  the  role  of  the  Hamiltonian  function.  If  the  flow  is 
time-periodic  then  the  study  of  (2)  is  typically  reduced  to  the 
study  of  a  two-dimensional  area  preserving  Poincare  map,  of 
which  the  standard  map  is  an  example.  Practically  speaking, 
the  reduction  to  a  Poincare  map  means  that  rather  than  viewing 
a  particle  trajectory  as  a  curve  in  continuous  time,  one  views 
the  trajectory  only  at  discrete  intervals  of  time,  where  the 
interval  of  time  is  the  period  of  the  velocity  field.  The  value  of 
making  this  analogy  with  Hamiltonian  dynamical  systems  lies 
in  the  fact  that  a  variety  of  techniques  in  this  area  have 
immediate  applications  to,  and  implications  for,  transport  and 
mixing  processes  in  fluid  mechanics.  For  example,  the  persist¬ 
ence  of  invariant  curves  in  the  Poincare  map  (KAM  curves) 
gives  rise  to  barriers  to  transport,  chaos  and  Smale  horseshoes 
provide  mechanisms  for  the  “randomization”  of  fluid  particle 
trajectories,  an  analytical  technique,  Melnikov’s  method,  al¬ 
lows  one  to  estimate  fluxes  as  well  as  describe  the  parameter 
regimes  where  chaotic  fluid  particle  motions  occur,  a  rela¬ 
tively  new  technique,  lobe  dynamics,  enables  one  to  efficiently 
compute  transport  between  qualitatively  different  flow  re¬ 
gimes,  bifurcation  theory  enables  one  to  understand  how 
qualitatively  different  flow  regimes  appear  and  disappear  as 
system  parameters  are  varied,  etc.  In  short,  many  new  methods 
for  the  study  of  fluid  transport  and  mixing  are  available  from 
dynamical  systems  theory.  In  the  following  we  will  show  how 
this  framework  and  techniques  can  be  used  in  a  problem  of 
geophysical  interest. 

The  techniques  and  results  described  in  the  following 
have  been  obtained  over  the  past  few  years  in  collaboration 
with  my  colleagues  Anthony  Leonard  and  John  Brady  at 
Caltech,  and  former  graduate  students  Vered  Rom-Kedar, 
Roberto  Camassa,  Tasso  Kaper,  Darin  Beigie,  and  Igor  Mezic', 

Transport  in  a  Cellular  Flow: 

The  Physical  Setting  and  the 
Model  Flow 

The  physical  setting  that  we  consider  is  as  follows.  We 
consider  a  (steady)  convection  cell  whose  horizontal  length  is 
much  larger  than  its  height  and  where  the  convection  cells  are 
aligned  along  the  y-axis  (Figure  3).  In  this  situation  the  flow 
is  essentially  two-dimensional  and,  assuming  stress-free 
boundary  conditions  and  single-mode  convection,  an  explicit 
form  of  the  velocity  field  is  given  by  (see  Chandrasekhar 
[1961]) 

i  =  -  cos(7rz)  sin(/::x)  =  {x,z),  (3) 

K  OZ 

z  =  Asin(7rz)cos(/^c)  =  ^{x,z),  (4) 
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where 


\l/(;c^)  =  j  sin{kx)  sliX^tz), 

A  is  the  maximum  vertical  velocity  in  the  flow, 


(K  the  wavelength  associated  with  the  cell  pattern),  and  the 
length  measures  have  been  non-dimensionaJized  so  that  the 
top  is  z=l  and  the  bottom  z=0.  This  flow  has  a  countable 
infinity  of  saddle  type  stagnation  points,  or  hyperbolic  fixed 
points  in  dynamical  systems  terminology,  on  the  upper  bound¬ 
ary  at 


(;c,z)  =  ('^,iy  =  0+l,±2,... 

and  a  countable  infinity  of  hyperbolic  fixed  points  on  the  lower 
boundary  at  (x,z)  =  (^,ov= 0^14:2.....  Fixed  points  with  the  same  x 

coordinate  are  connected  by  a  stagnation  streamline,  or  het¬ 
eroclinic  orbit  in  dynamical  systems  terminology.  The  result 
is  an  infinite  number  of  cells,  and  many  of  the  fluid  mechanical 
problems  of  interest  are  concerned  with  the  transport  of  a 
passive  tracer  from  cell-to-cell.  Fluid  flows  which  are  divided 
up  into  a  collection  of  cells  occur  in  a  variety  of  geophysical 


Figure  2. 

Stable  and  unstable  manifolds  of  a  period  one  point  in  a  chaotic 
region. 


applications  and  the  methods  and  approaches  discussed  in  this 
article  in  principle  can  be  applied  to  any  such  cellular  fiow. 

The  Experiments  of  Solomon 
and  Gollub 

We  motivate  the  analysis  by  discussing  some  experiments 
of  Solomon  and  Gollub  on  this  system.  If  the  temperature 
difference  between  the  lop  and  bottom  of  the  convection  cell 
is  increased  an  additional  time-periodic  instability  occurs, 
resulting  in  a  time-periodic  velocity  field  (Clever  and  Busse 
[1974],  Bolton  et  ai  [1986]).  Solomon  and  Gollub  [1988] 
experimentally  studied  the  transport  of  dye  in  such  a  flow  and 
they  made  the  following  observations. 

1.  There  was  a  dramatic  enhancement  of  the  “effective 
diffusivity”. 

2.  Molecular  diffusion  appeared  to  play  no  role  in  the 
transport. 

3.  The  flux  across  the  cell  boundaries  (heteroclinic  orbits) 
depended  linearly  on  the  amplitude  of  the  oscillatory 
instability  and  was  independent  of  the  wavelength  X. 

Therefore,  the  transport  problem  is  radically  different  in 
the  unsteady  case  as  compared  to  the  steady  case  since  in  the 
steady  case  streamlines  cannot  cross  and  therefore  tracer 
crosses  the  cell  boundaries  solely  as  a  result  of  molecular 
diffusion.  In  order  to  study  these  issues  Solomon  and  Gollub 
introduced  the  following  model  of  the  even  oscillatory  cell 
instability: 

X  =  -^cos{Kz)[sin{kx)  +  zkf{t)cQskx)]  =  -^(XyZ,t) 

z  =  Asin(7rz)[cos(/::x)  --  Elfit)sir\{kx)]  =  (5) 

where 


\\f(x,z,t)  ~—sin(kx)sm(nz)  +  8y(0cos(^x)sin(7tz) 


and  we  will  take 

fit)  =  cos(co^)  and  e  (R  -  R^)  , 

where  R  is  the  Rayleigh  number  and  R^  its  critical  value  at 
which  the  time-periodic  instability  occurs.  The  pros  and  cons 
of  this  model  are  discussed  in  Solomon  and  Gollub  [1988]. 
This  model  will  be  the  starting  point  of  our  analysis  and, 
motivated  by  the  experimental  observations  of  Solomon  and 
Gollub,  we  will  consider  the  following  four  questions, 

1.  What  is  the  mechanism  for  cell-to-cell  transport? 

2.  Can  we  quantify  the  spreading  of  a  passive  tracer  (dye) 
initially  contained  in  one  cell? 
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Figure  3. 

Streamlines  for  the  steady,  cellular  flow. 


3.  Can  we  predict  the  number  of  cells  invaded  as  a  function 
of  time? 

4.  What  are  the  effects  of  the  addition  of  a  small  amount 
of  molecular  diffusion? 

The  Mechanism  for  Cell-to-Cell 
Transport 

In  the  case  of  unsteady  two-dimensional  flows  modem 
dynamical  systems  theory  provides  us  with  some  new  methods 
and  concepts  that  enable  us  to  analyze  the  particle  kinematics. 
For  example,  the  saddle-type  stagnation  points  may  move  and 
the  separatrices  associated  with  different  “moving  stagnation 
points”  typically  will  not  coincide  to  form  cell  boundaries  in 
the  flow.  This  situation  can  lead  to  cell-to-cell  transport.  Very 
roughly  speaking,  dynamical  systems  theory  provides  us  with 
methods  for  studying  particle  kinematics  in  unsteady  flows 
with  moving  stagnation  points  and  their  associated  moving 
separatrices.  For  the  case  of  time-periodic  variability,  if  the 
variability  is  weak  then  it  can  be  shown  that  a  saddle-type 


stagnation  point  persists,  but  it  oscillates  periodically  in  time 
(with  the  period  of  the  variability).  Moreover,  this  “oscillating 
stagnation  point”  (more  precisely,  hyperbolic  periodic  orbit) 
still  maintains  its  saddle-type  stability  characteristic  and  asso¬ 
ciated  with  it  are  two  one-dimensional  periodically  varying 
curves  called  the  stable  and  unstable  manifolds  of  the  hyper¬ 
bolic  periodic  orbit.  Trajectories  on  the  stable  manifold  ap¬ 
proach  the  hyperbolic  periodic  orbit  asymptotically  as  r  -> 
and  trajectories  on  the  unstable  manifold  approach  the  hyper¬ 
bolic  periodic  orbit  asymptotically  as  r  -  oo.  These  curves 
are  invariant  curves  in  the  sense  that  particle  trajectories  that 
start  on  the  curves  remain  on  the  curves  for  all  time.  This 
property  of  invariance  implies  that  other  trajectories  cannot 
cross  these  curves,  hence  they  act  as  barriers  to  fluid  transport 
in  a  way  that  we  will  describe  more  fully  shortly.  These  curves 
are  the  unsteady  analogs  of  the  separatrices  that  commonly 
arise  in  steady  flows.  However,  as  opposed  to  the  steady  case, 
the  stable  and  unstable  manifolds  of  different  oscillating  sad¬ 
dle-type  stagnation  points  typically  do  not  coincide  to  form 
complete  barriers  to  fluid  transport.  Rather,  at  a  fixed  instant 
of  time  they  may  intersect  in  a  discrete  set  of  points.  Never¬ 
theless,  these  curves  still  have  the  important  kinematical  effect 
of  forming  the  boundaries  between  fluid  particle  trajectories 
that  move  in  different  directions  after  their  encounter  with  the 
oscillating  saddle-type  stagnation  point.  However,  since  the 
stable  and  unstable  manifolds  of  different  oscillating  saddle- 
type  stagnation  points  need  not  join  up  (and,  hence,  form  a 
curve  of  finite  length  in  the  flow),  they  may  have  infinitely 
long  length  and  wind  throughout  the  flow  (and  themselves)  in 
a  complicated  geometrical  fashion.  Keeping  in  mind  the  kine¬ 
matical  meaning  of  these  curves  as  boundaries  between  quali¬ 
tatively  different  fluid  particle  motions,  they  therefore  act  as 
partial  barriers  to  fluid  transport.  As  a  result,  fluid  transport 
in  two-dimensional  unsteady  flows  may  be  much  more  varied 
and  complex.  Dynamical  systems  theory  gives  us  a  way  to 
make  these  intuitive  notions  precise  and  quantitative,  which 
we  now  describe  for  time-periodic  variability. 

For  time-periodic  variability  a  variety  of  standard  dy¬ 
namical  systems  techniques  and  results  make  it  advantageous 
to  consider  particle  kinematics  via  the  so  called  Poincare  map, 
which  we  will  denote  by  /.  That  is,  rather  than  plotting  a 
particle  trajectory  as  a  continuous  curve  in  space  one  plots  its 
evolution  at  discrete  intervals  of  time,  where  the  interval  of 
time  is  equal  to  the  period  of  the  variability.  In  this  setting,  the 
oscillating  saddle-type  stagnation  point  is  manifested  as  fixed 
point  of  the  Poincare  map  and  the  stable  and  unstable  mani¬ 
folds  associated  with  the  fixed  point  are  similarly  fixed  in 
space  for  the  Poincare  map.  Consequently,  the  situation  bears 
some  similarity  to  the  steady  case,  but  it  is  important  to  keep 
in  mind  that  for  the  Poincare  map  particle  trajectories  are 
manifested  as  sequences  of  discrete  points,  rather  than  con¬ 
tinuous  curves.  Also,  as  we  mentioned  earlier,  the  stable  and 
unstable  manifolds  of  different  hyperbolic  fixed  points  may 
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Figure  4. 

The  behavior  of  the  stable  and  unstable  manifolds  on  the  walls  of 
the  convection  cell  under  the  time-periodic  perturbation-  the 
heteroclinic  tangle,  numerically  calculated  forz  =  0.1,  (a  =  0.6, 
A  =  0.1. 
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not  coincide(or  “join  up”)  to  make  a  separatrix  of  finite  length, 
rather  they  may  intersect  in  a  discrete  number  of  points  and 
wind  throughout  the  flow.  We  demonstrate  this  for  (5)  in 
Figure  4.  The  stable  and  unstable  manifolds  of  the  different 
hyperbolic  fixed  points  actually  have  infinite  length,  but  we 
obviously  can  only  show  a  finite  portion  of  each  in  the  figure. 
Next  we  show  how  the  geometric  structure  associated  with  the 
stable  and  unstable  manifolds  of  the  different  hyperbolic  fixed 
points  influence  transport. 

For  the  Poincar^  map,  segments  of  finite  length  of  the 
stable  and  unstable  manifolds  of  the  hyperbolic  fixed  points 
can  be  used  to  form  cell  boundaries.  We  illustrate  this  in  Figure 
5  where  the  segments  of  finite  length  begin  at  the  hyperbolic 
fixed  points  and  end  at  particular  intersection  points  of  the 
stable  and  unstable  manifolds,  which  we  will  refer  to  as  the 
boundary  intersection  points.  Now  we  can  consider  the  ques¬ 
tion  of  fluid  exchange  between  the  different  regions  defined 
by  our  chosen  boundaries  that  are  shown  in  Figure  5. 

First  we  list  two  rules  that  points  on  the  stable  and  unstable 
manifolds  of  hyperbolic  fixed  points  must  obey  under  iteration 
by  the  Poincare  map.  (By  “iteration”  we  mean  repeated  appli¬ 
cation  of  the  Poincare  map  to  a  particle  or,  in  other  words, 
discrete  time  evolution  of  a  particle  trajectory.) 

1 .  A  point  that  is  on  both  the  stable  and  unstable  manifolds 
must  remain  on  both  manifolds  under  all  positive  and 
negative  iterations  of  the  Poincare  map.  This  is  because 
these  manifolds  are  invariant  manifolds. 

2.  Points  on  the  stable  or  unstable  manifolds  of  a  fixed 
point  maintain  their  relative  order  (in  the  sense  of  dis¬ 
tance  in  arc  length  from  the  fixed  point)  under  iteration 
by  the  Poincare  map.  That  is,  if  we  consider  two  distinct 
points  on  the  stable  or  unstable  manifold  then  one  will 
be  closer  than  the  other  to  the  fixed  point  in  the  sense  of 
arc  length  along  the  stable  or  unstable  manifold.  If  we 


then  iterate  both  points,  the  same  point  will  always 
thereafter  be  closer  to  the  fixed  point.  The  reasons  for 
this  are  more  technical  and  are  related  to  uniqueness  of 
particle  trajectories  for  a  given  initial  condition. 

With  these  two  rules  in  mind  let  us  return  to  Figure  5  and 
consider  the  preimages  under  the  Poincare  map  of  the  bound¬ 
ary  intersection  points,  i.  e. ,  the  backwards  time  evolution  of 
these  points  for  one  period  of  the  variability.  By  rule  1  these 
points  still  lie  on  both  the  stable  and  unstable  manifolds,  and 
since  the  manifolds  are  smooth  curves  they  must  wind  through 
each  other  as  shown  in  Figure  6.  In  this  figure  we  have  shown 
also  precisely  one  additional  point  of  intersection  of  the  mani¬ 
folds  between  the  boundary  intersection  point  and  its  pre¬ 
image.  There  is  a  technical  reason  for  this  that  we  briefly  must 
mention.  Generally,  there  must  be  an  odd  number  of  additional 
intersection  points  between  the  boundary  intersection  point 
and  its  preimage  if  we  have  uniqueness  of  parcel  trajectories 
for  a  given  initial  condition.  Without  loss  of  generality,  we  can 
assume  that  we  have  one  additional  intersection  point  since  all 
other  cases  can  be  recast  in  this  form.  Hence,  the  segments  of 
the  stable  and  unstable  manifolds  between  a  boundary  inter¬ 
section  point  and  its  preimage  form  two  lobes  and  these  two 
lobes  are  said  to  form  a  turnstile.  The  turnstiles  are  the  mecha¬ 
nisms  governing  transport  between  the  different  regions,  as  we 
will  now  argue.  Consider  the  image  of  the  turnstile  lobes  under 
the  Poincare  map  f  (i.  e. ,  the  evolution  of  these  points  for  one 
period  of  the  flow).  Using  the  two  rules  above,  as  well  as  the 
fact  that(for  a  continuous  and  continuously  invertible  map)  the 
boundary  of  a  set  maps  to  the  boundary  of  its  image  and  the 
interior  of  a  set  maps  to  the  interior  of  its  image,  we  see  that 
the  image  of  the  turnstile  lobes  appear  as  in  Figure  6.  Thus, 
the  respective  turnstile  lobes  have  “switched  sides”  in  relation 
to  the  boundary  (hence  the  term  “turnstile”).  Now,  using  rule 
2  above,  it  can  be  argued  that  the  only  points  that  cross  the 
boundary  in  one  iteration  of  the  Poincare  map  are  those  in  the 


Figure  5. 

Cell  boundaries  for  the  unsteady  flow(for  the  associated  Poin¬ 
care  map). 
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Figure  6. 

The  turnstile  lobes  associated  with  cell-to-cell  transport. 
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turnstile  lobes.  Therefore  the  flux  across  the  boundary  in  one 
period  is  simply  the  area  of  the  turnstile  lobe. 

In  the  case  where  the  variability  is  weak  there  is  a  pertur¬ 
bation  method  that  will  enable  us  to  estimate  the  lobe  area-and 
therefore  the  flux  per  period-across  the  “broken  separatrices”. 
This  is  Melnikov’s  method  for  time-periodic  vector  fields.  The 
Melnikov  function  is  a  measure  of  the  distance  between  the 
stable  and  unstable  manifolds  measured  along  a  line  perpen¬ 
dicular  to  the  unperturbed  separatrix  at  the  point  parametrized 
by  to,  where  tg  can  be  thought  of  as  an  arc  length  parameter 
measuring  distance  along  the  unperturbed  separatrix.  In 
Camassa  and  Wiggins  [1991]  the  Melnikov  function  on  the 
zero  phase  cross  section  is  calculated  and  found  to  be  given 
by: 

M{t^  =  (x^sech{-^)sm{mQ), 

Clearly  the  Melnikov  function  has  a  countable  infinity  of 
simple  zeros,  that  correspond  to  a  countable  infinity  of  inter¬ 
sections  of  the  stable  and  unstable  manifolds. 

As  described  above,  the  flux  of  fluid  between  cells  is 
given  by  the  area  of  the  turnstile  lobes.  In  this  perturbative 
setting  an  approximation  to  the  flux  (lobe  area)  is  given  by  the 
integral  of  the  Melnikov  function  between  two  adjacent  zeros, 
i.  e. , 

(i) 

M-(^i,o(l))  =  4(A).i(1))  =  2ejec/^  +  cs{z\ 

where  (L^  j(l))  denotes  the  area  of  the  lobe  L;  j.  We  remark 
that  because  of  translational  symmetries,  this  result  holds  for 
all  turnstile  lobes.  The  symmetries  of  this  problem  are  dis¬ 
cussed  in  Camassa  and  Wiggins  [1991].  This  approximation 
to  the  flux  show  us  that  the  flux  depends  linearly  on  the 
amplitude  of  the  oscillatory  instability  and  is  independent  of 
the  wavelength  of  the  cell  patterns  X,  exactly  as  observed 
experimentally  by  Solomon  and  Gollub. 


The  Number  of  Cells  Invaded 
as  a  Function  of  Time:  Lobe 
Dynamics 

We  next  provide  another  example  of  how  “lower  dimen¬ 
sional  structures”  in  the  flow,  the  lobes,  determine  a  large  scale 
transport  property  of  the  flow.  Suppose  a  given  cell  is  filled 
with  a  passive  tracer.  How  many  cells  are  “invaded”  by  this 
tracer  after  some  elapsed  time?  The  answer  to  this  question 
can  be  determined  by  the  geometry  of  the  lobe  intersections 
associated  with  the  cell  boundaries.  Using  lobe  dynamics  type 
arguments  and  invariance  of  the  manifolds,  as  well  as  transla¬ 
tion  and  reflection  symmetries,  an  upper  and  lower  bound  can 
be  placed  on  the  time  that  a  passive  scalar  initially  in  Rj  first 
invades  R.j.  Surprisingly,  all  the  necessary  information  is 
contained  in  two  integers  which  we  refer  to  as  the  signatures 
of  the  heteroclinic  tangles  forming  the  cell  boundaries, 
namely: 

•  m  =  the  smallest  integer  such  that/"(Lj  Q(i))nLQ_j(i)4). 

•  m'  =  the  smallest  integer  such  that  the  boundary  of 

intersects  the  boundary  ofL.j  ,2(1)  in  four 
distinct  points 

The  signatures  are  illustrated  in  Figure  7  form  =  1  and  m'  =  2. 

If  we  denote  the  period  of  the  velocity  field  by 
7^—  and/j^  the  time  necessary  for  tracer  initially  in  Rj  to  enter 
R_j  we  have  the  following  general  result  obtained  in  Camassa 
and  Wiggins  [1991]: 

(jrh+\)T<  tf  <  [(j -  l)rn  +(m+  l)]r,  (j >  2). 

Therefore,  the  upper  and  lower  bounds  for  the  first  invasion 
time  are  completely  determined  by  the  signatures  m  and  m'.  This 
result  is  independent  of  the  boundary  conditions.  For  instance. 


Figure  7. 

An  illustration  of  the  signatures  of  the  heteroclinictangle  asso¬ 
ciated  with  the  cell  boundaries. 
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numerical  computations  show  that  in  the  case  A=0. 1,  e  =  0. 1 
and  (0  =0.6  we  have  mand^=3  indicating  that  one  cell  is 
invaded  every  three  periods,  for  A=0.1,  e  =  0.1  and  co  =  0.24 
we  have  mand^  =1  indicating  that  one  cell  is  invaded  every 
period,  and  for  A=0. 1 ,  z  =0.01  and  co=0. 6  we  have  m  and  m'.  = 
4  indicating  that  one  cell  is  invaded  every  four  periods.  We 
remark  that  the  signatures  are  very  easy  to  compute  numeri- 
cally  as  one  needs  to  compute  only  a(usually  short)  finite 
length  of  the  manifolds.  Moreover,  the  necessary  manifold 
interections  are  typically  robust  with  respect  to  numerical 
computations. 

Quantifying  the  Spreading  of  a 
Passive  Scalar:  Lobe  Dynamics 

The  dynamics  of  the  turnstile  lobes,  lobe  dynamics,  can 
be  used  to  rigorously  determine  the  amount  of  tracer  originally 
in  one  cell  that  is  contained  in  any  other  cell  at  a  later  time. 
Suppose  the  cell  Rj  is  uniformly  filled  with  tracer  (species  S^) 
at  t=0  how  much  of  species  Sj  is  contained  in  R-  at  any  t  >  0? 
If  we  denote  by  j  (n)  the  total  amount  of  species  Sj  in  Rj 
immediately  after  the  n-th  iterate  and  by  a^^  j(n)  =Ti  j(n)-Ti 
j(n-l)  the  flux  of  species  into  Rj  on  the  n-*th  iterate  ,  then  it 
can  be  shown  that  this  last  formula  can  be  reduced  to  an 
expression  containing  areas  of  intersection  sets  involving  im¬ 
ages  of  only  one  of  the  turnstile  lobes  q(1)  associated  with 
the  boundary  of  Rj.  The  resulting  expression  is  given  by 

aij(n)  =  Tijin)  -  TJ(n  -  1)  =  (5, 2  +  djPML^  oH)) 

n-1 

+ 1  o(l))  j 

where,  recall,  |i  (L^  j(l))  denotes  the  area  of  the  lobe  j(l). 
Consequently,  a  relatively  straightforward  numerical  compu¬ 
tation  of  j(n)  is  possible.  To  compute  a^  j(n)  one  merely 
locates  the  lobes  indicated  in  (6),  iterates  Lj  q(1)  (i.  e. ,  a  grid 
of  points  covering  Lj  ^(l)),  determines  the  area  of  intersection 
of  the  appropriate  iterates  of  Lj  o(l)  with  the  other  lobes  in  (6), 
and  adds  up  the  result  according  to  the  recipe  provided  by  (6). 
Numerically  this  method  is  extremely  efficient  and  affords 
tremendous  savings  in  CPU  time  over  standard  “brute  force” 
methods;  details  are  discussed  in  Camassa  and  Wiggins 


[1991].  However,  practically,  this  computation  can  only  be 
performed  for  a  finite  amount  of  time.  In  the  next  section  we 
discuss  an  infinite  time  result. 

Dispersion 

Dispersion  is  a  quantity  that  is  frequently  of  interest  in 
geophysical  settings.  For  example,  Pasmanter  [1991]  and  Rid- 
derinkhof  and  Zimmerman  [1992]  have  recently  studied  dis¬ 
persion  in  models  of  shallow  tidal  flows  and  found  that  it 
behaves  asymptotically  like  Weiss  and  Knobloch  [1989] 
have  performed  a  numerical  study  of  dispersion  in  modulated 
traveling  waves  in  binary  fluid  convection  and  found  a  disper¬ 
sion  exponent  of  1.93.  As  their  results  were  numerical,  they 
were  only  able  to  compute  for  a  finite  length  of  time.  Using  a 


Figure  8. 


Comparison  of  lobe  dynamic  transport  results  (solid)  with  trans¬ 
porttaking  into  effect  molecular  diffusivity  (dashed)  for  the  three 
cases.  a)A=  0.1,  o>  =  0.6,  z  ~  0.1,  b)A=  0.1,  (o  =  0.6,  e  =  0.01. 
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simple  application  of  Birkhoff ’s  ergodic  theorem  we  can  de¬ 
rive  a  rigorous  result  that  applies  in  each  of  these  settings  and 
gives  conditions  under  which  the  dispersion  should  behave 
asymptotically  (in  time)  like  As  such,  it  may  also  prove  to 
be  a  useful  guide  for  numerical  investigations.  Our  discussion 
will  be  in  the  context  of  the  time-periodic  cellular  flow,  how¬ 
ever  the  results  can  be  applied  more  generally  (See  Mezic'  and 
Wiggins  [1994b],  [1994c]). 

Referring  to  the  velocity  field  (5),  the  displacement  of  a 
particle  trajectory  in  the  x  direction  is  given  by 

jc(0  -x(0)  =  f  (-^cos(nz{t))[sm{kx{t))  + 

(7) 

tk  cos(co0cos(fcc(0)])rf^, 

where  (x(t),  z(t))  denotes  a  trajectory.  For  the  application  of 
Birkhoff’s  ergodic  theorem  we  need  to  introduce  a  standard 
dynamical  systems  “trick”  that  is  applied  to  time-periodic 
vector  fields.  That  is,  we  introduce  the  phase,  call  it  0,  of  the 
time-periodic  part  as  a  new  dependent  variable  and  append  it 
to  the  vector  field.  The  new  component  is  therefore  6  =  od  and 
the  vector  field  is  then  viewed  as  a  three-dimensional  time-in- 
dependent  dynamical  system.  If  D  denotes  the  flow  domain 
(periodic  in  x,  bounded  in  z),  then  the  phase  space  of  this 
“enlarged”  dynamical  system  is  denoted  D  x  S\  where 
denotes  the  phase  of  the  time  periodic  part.  The  “flow”  gener¬ 
ated  by  this  dynamical  system  is  then  denoted  by  0^(x,z,9),  i. 
e. ,  this  is  the  trajectory  that  starts  at  the  point  (x,  z,  0).  If  we 
denote  the  x  component  of  the  velocity  of  (5)  by  m,  then  (7) 
can  be  succinctly  written  as 

;c(0  -  x(0)  =  /^«((p/x,z,e))d;,  (8) 


We  are  interested  in  determining  the  asymptotic  behavior 
of  the  dispersion.  We  have  the  following  calculations: 


t->oo  t 

^  ^  \  2 
=  lim  /  (-  }  M  ((p.j(x,z,0))  c/t  -  <  I  j  u{%{x^ff))d\ ) )  ), 

'o 

=  <  (lim  -  J  u  ((p^(x,z,e))d'i  -  <  lim-  J  u((p^(.x^fi))dx) )  >. 

=  ((uV^,0)- 

1 .  The  passage  from  the  first  to  the  second  line  is  justified 
by  the  fact  that  the  function  u  is  bounded  and  integrable 
on  D  X  In  the  second  line  the  limit 

2.  In  the  second  line  the  limit 

t 

lim- jt/((jV(x,z,0))tir  = 

exists  for  all  points  in  D  x  by  Birkhoff ’s  ergodic 
theorem  (see  Arnold  and  Avez  [1968]),  with  the  possible 
exception  of  a  set  of  \x-measure  zero.  This  limit  is  the 
time  average  of  the  function  w  along  the  fluid  particle 
trajectory  that  starts  at  the  point  (x,  z,  0).  Moreover, 
Birkhoff’s  ergodic  theorem  also  guarantees  that  this 
limit  is  integrable.  This,  together  with  the  boundedness 
of  M,  implies  that  the  quantity  a  defined  above  is  finite, 
and  if  it  is  nonzero,  we  can  conclude  that  the  dispersion 
of  the  ensemble  of  particles  in  the  z  direction  behaves 
asymptotically  like  t^. 


Figure  9. 

Contours  corresponding  to  points  having  the  same  timeaver- 
age  along  trajectories  for  the  x-component  of  velocity,  u . 


The  mean  square  displacement  or  dispersion  of  the  x 
component  of  (5)  of  an  ensemble  of  points  under  the  flow  is 
given  by 


((x(t)  -  x(0)  -  (x(t)  -  x(0)))^)  =  D^(t), 

where  the  average  indicated  by  the  angle  brackets  is  defined 
as 


((x(0-x(0))^  i  ^{x{t)-x{0))pdii, 

DxS^ 

p  =  p(x,  z,  0)  is  the  initial  distribution  of  points  (assumed  to  be 
bounded  and  integrable  on  £>  x  S  ^),  and  dp,  denotes  the  measure 
or  “volume  element”  on  D  x  S*'  Incompressibility  of  the  flow 
implies  that  the  flow  is  “measure  preserving”,  which  is  impor¬ 
tant  for  the  application  of  Birkhoff’s  ergodic  theorem. 
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Figure  10. 
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The  nature  of  the  coefficient  a  gives  some  insight  into  the 
dynamical  mechanism  giving  rise  to  t^  dispersion.  It  is  easy  to 
see  that  since  the  expression  inside  the  angle  brackets  defining 
a  is  nonnegative,  a=0  if  and  only  if  u*  =  <u*>  on  the  support 
of  /?,  i.  e.,  on  the  set  of  points  for  which  p(x,  z,  9)  is  nonzero, 
with  the  possible  exclusion  of  sets  of  measure  zero.  However, 
<u*>  is  a  constant.  Therefore,  we  can  make  the  following 
conclusion.  Let  C  a  denote  the  support  of  p(x,  z,  0). 
Then,  if  m*(x,  z,  0)  is  not  constant  almost  everywhere  on  C, 
t^  as  t  — >  oo.  Now  assume  that  the  flow  is  ergodic.  Then 
u*=  <u*>  =  <u>  almost  everywhere.  Therefore,  a=0.  So,  a 
necessary  condition  for  t^  dispersion  is  non-ergodicity  of  the 
flow. 

We  end  with  some  final  observations. 

1.  Note  that  our  result  is  independent  of  the  Reynolds 
number.  We  are  dealing  solely  with  kinematical  consid¬ 
erations. 

2.  There  has  been  much  work  done  in  the  last  10  years  on 
chaotic  fluid  particle  dynamics  in  two  dimensional, 
time-periodic  velocity  fields.  In  such  situations  one 
typically  sees  a  mixture  of  regular  and  chaotic  regions 
of  fluid  particle  motions  in  the  flow.  Our  dispersion 
result  implies  that  if  one  takes  an  initial  distribution  of 
points  that  is,  roughly  speaking,  not  entirely  contained 
in  a  single  regular  or  chaotic  region,  then  the  asymptotic 
behavior  of  the  dispersion  will  go  like  t^. 

3.  The  case  of  a  diffusive  tracer  has  recently  been  treated 
by  Mezic'  et  al.  [1994].  This  work  shows  that  the  t^ 
dispersion  result  in  the  non-diffusive  case  is  important 
for  studying  the  diffusive  case  in  the  high  Peclet  number 
limit.  (The  Peclet  number  is  a  dimensionless  number 
that  is  typically  of  the  form  Pe=^,  where  U  is  a  charac¬ 


teristic  velocity,  I  is  a  characteristic  length,  and  D  is  the 
diffusivity.)  In  particular,  for  the  time-periodic  cellular 
flow  it  is  possible  to  show  that  t^  dispersion  of  the 
convective  part  of  the  problem  implies  a  Pe^  dependence 
of  the  effective  diffusivity  coefficient  in  the  diffusive 
case.  In  fact,  the  dependence  of  the  effective  diffusivity 
on  the  Peclet  number  can  be  determined  solely  on  the 
knowledge  that  the  convective  part  of  the  problem  ex¬ 
hibits  t^  dispersion.  Numerical  studies  of  chaotic  incom¬ 
pressible  flows  typically  show  regions  of  chaotic  and 
ordered  behavior,  therefore  non-ergodicity.  In  order  to 
study  dispersion  in  chaotic  regions,  the  usual  procedure 
is  to  consider  an  initial  distribution  of  points  that  is 
entirely  contained  in  what  seems  to  be  a  chaotic  (and 
supposedly  ergodic)  region.  From  the  above  result,  in 
the  diffusive  case  it  would  not  make  any  difference 
whether  the  points  were  initially  placed  in  only  one 
ergodic  region,  or  spread  out  over  several  ergodic  re¬ 
gions.  The  only  importance  lies  in  the  fact  that  in  the 
related  non-diffusive  problem,  when  the  particles  are 
placed  initially  in  several  ergodic  regions,  the  nondif- 
fusive  dispersion  behaves  like  t^  at  large  times.  Then  the 
Pe^  regime  is  obtained  in  the  diffusive  problem.  In  this 
context  we  consider  these  results  to  be  typical  fora  large 
class  of  laminar  flows. 

Relative  Time  Scales  of  Lobe 
Transport  Versus  TVansport  by 
lecular  Diffusion  through  Lobe 
Geometiy 

So  far,  we  have  seen  how,  in  the  absence  of  molecular 
diffusion,  lobe  dynamics  can  be  used  to  quantify  the  spreading 
of  a  passive  tracer  for  finite  times.  Also,  we  have  seen  how 
ergodic  theory  can  enable  us  to  understand  the  asymptotic  (in 
time)  behavior  of  dispersion.  In  this  section  we  show  how  the 
lobes  can  be  used  to  determine  a  time  scale  over  which  lobe 
transport  dominates  molecular  diffusion.  Our  knowledge  of 
the  dynamics  of  fluid  particles  associated  with  lobes  suggests 
the  following  criterion  that  might  explain  the  time  scale  over 
which  lobe  transport  dominates  molecular  diffusion. 

The  time  scale,  T^,for  a  tracer  to  diffuse  across  a  distance  of 
the  order  of  a  turnstile  width  should  be  long  compared  to  the 
time  it  would  take  for  a  lobe  to  be  mapped  across  the  boundary 
of  a  region,  i.  e.  T. 

Hence,  we  have 

r 

where  3(8)  is  the  maximum  width  of  a  turnstile  lobe.  Using  the 
a(e)  approximation  to  the  turnstile  width  given  by  the  Mel¬ 
nikov  function,  we  have 
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According  to  our  criterion,  lobe  transport  will  dominate 
diffusive  transport  provided  T^j »  T.^We  check  this  for  certain 
parameter  values.  With  z)  =  5.0  x  and  A=  0. 1  we  have  the 
following  three  cases: 

1 .  03  =  0.6, e  =  0.1  implies  To  *  2007, 

2.  03  =  0.24,e  =  0.1  implies  To  ~  3007, 

3.  03  =  0.6,e  =  0.01  implies  To  «  27, 

Thus,  in  the  first  two  cases  we  would  expect  lobe  transport 
to  dominate  molecular  diffusion  for  about  200  and  300  peri¬ 
ods,  respectively,  and  in  the  last  case  only  for  about  2  periods. 
This  can  be  check  by  considering  the  generalized  Langevin 
equations,  which  describe  the  influence  of  molecular  diffusiv- 
ity: 


i=-^+n«). 

(9) 

(10) 

where  il(0>C(0  (diffusion  terms)  are  random  variables  with  a 
Gaussian  probability  distribution  ,  such  that  their  correlations 
are: 

<n(0n(O  > = <i;(oC(^0>=2D5(z  - (i  i) 
<n(OC(0>=o.  (12) 

In  Figure  8  we  compare  the  results  of  the  lobe  dynamics 
calculations  (solid  lines),  which  were  obtained  by  numerically 
implementing  the  procedure  described  in  Section  5,  and  the 
calculations  which  consider  molecular  diffusion(dashed 
lines),  which  were  obtained  by  integrating  (10).  From  this 
figure  we  can  see  that  our  criterion  works  reasonably  well. 
Thus  we  see  that  geometrical  and  dynamical  features  of  the 
lobes  have  lead  us  to  a  physical  argument  that  gives  us  a  time 
scale  on  which  molecular  diffusion  can  be  neglected  in  the 
process  of  inter-cell  transport. 

“Patchiness”  in  Fluid  Flows 
and  Distribution  Functions 

Finally,  we  want  to  describe  a  new  phenomenon  that  may 
prove  to  be  a  useful  way  of  describing  flows  in  a  variety  of 
geophysical  settings.  Pasmanter  [1988],  [1991]  has  studied 
mechanisms  that  give  rise  to  the  variability  of  dispersion 
processes  in  the  ocean.  A  particularly  important  phenomenon 
to  which  he  referred  is  known  as  “patchiness”,  i.  e. ,  a  situation 


where  parts  of  a  distribution  of  passive  tracer  may  disperse  at 
different  speeds  compared  to  its  surroundings.  We  want  to 
show  that  the  mathematical  framework  developed  in  Mezic' 
and  Wiggins  [1994c]  can  be  useful  for  studying  this  phenome¬ 
non.  We  will  illustrate  this  with  the  cellular  flow  example. 

In  the  cell-to-cell  transport  problem  the  time  average 
along  particle  trajectories  of  the  x-component  of  velocity  is  a 
useful  quantity.  In  order  to  illustrate  the  “patchiness”  effect  we 
numerically  determine  the  distribution  of  average  x-compo- 
nents  of  velocity  for  fluid  particle  trajectories  in  one  cell.  In 
Figure  9  we  plot  contours  corresponding  to  initial  points  of 
trajectories  that  have  the  same  average  x-velocity  which  we 
obtained  by  numerically  computing  these  time  averages  for  a 
uniform  grid  of  1600  points  for  a  time  interval  of  5000  periods. 
Regions  of  nonzero  average  x  velocity  correspond  to  “accel¬ 
erator  modes”,  i.  e.,  regions  of  initial  conditions  corresponding 
to  trajectories  with  nonzero  average  velocity  in  the  x  direction, 
and  these  are  the  points  that  participate  in  cell-to-cell  transport. 
In  Figure  10  we  plot  a  histogram  of  this  data,  i.  e. ,  we  plot  the 
number  of  points  corresponding  to  a  given  “bin  width”  of 
average  x  velocity,  where  the  bin  width  is  0.  002.  “Patches” 
are  mathematically  defined  as  sets  of  positive  measure  in  the 
flow  on  which  the  x  component  of  velocity  has  a  constant 
time-average  along  fluid  particle  trajectories.  In  Mezic'  and 
Wiggins  [1994c]  a  distribution  function  is  derived  that  allows 
one  to  deduce  certain  properties  related  to  the  patchiness.  In 
particular,  the  “patches”  correspond  to  points  of  discontinuity 
of  a  certain  distribution  function.  An  additional  result  obtained 
through  study  of  this  distribution  function  tells  us  that  there 
are  either  a  finite  number  of  patches,  or  there  is  a  countable 
number  of  them,  but  also  countably  many  of  them  are  of 
smaller  measure  than  any  prescribed  number.  In  Figure  9,  the 
“patches”  are  represented  by  the  darkest  and  brightest  spots. 
We  could  also  define  “patchiness”  as  the  regime  in  which  more 
than  one  “patch”  exists.  It  then  turns  out  that  non-ergodicity 
of  the  flow  is  a  necessary  condition  for  “patchiness”.  More 
precisely,  it  can  be  shown  that  if  the  flow  is  ergodic,  there  is 
only  one  “patch”,  represented  as  one  and  only  point  of  discon¬ 
tinuity  of  the  distribution  function  derived  in  Mezic'  and 
Wiggins  [1994c].  Figure  10  represents  the  “generalized  de¬ 
rivative”  of  this  distribution  function,  that  is,  the  peaks  corre¬ 
spond  to  different  “patches”.  From  these  results  one  can 
deduce  a  criteria  for  determining  when  a  flow  will  be  “patchy” 
If  the  flow  is  such  that  it  has  an  ergodic  component  of  positive 
measure,  and  if  the  initial  distribution  of  passive  tracer  is  such 
that  the  measure  of  that  component  is  also  of  positive  measure, 
then  the  flow  will  evolve  into  a  “patchy”  flow. 


Summary 

In  this  article  we  have  taken  a  specific  flow  and  shown 
how  an  arsenal  of  newly  developed  techniques  from  dynami¬ 
cal  systems  theory  can  be  brought  to  bear  on  the  study  of  the 
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Lagrangian  transport  properties  of  this  flow.  This  dynamical 
systems  approach  shows  much  promise  in  dealing  with  the 
types  of  Lagrangian  transport  issues  that  arise  in  geophysical 
fluid  dynamics.  Much  of  the  promise  lies  in  the  fact  that  the 
mathematical  framework  of  dynamical  systems  theory  is  re¬ 
markably  similar  to  the  experimental  framework  of  modem 
flow  visualization  (from,  e.  g. ,  tracer  data  or  satellite  obser¬ 
vations)  in  that  each  is  concerned  with  the  motion  in  time  over 
large  regions  of  space  and  the  role  that  lower  dimensional 
organized  “stmctures”,  or  geometrical  features,  in  space  play 
in  governing  the  motion.  Dynamical  systems  theory  provides 
a  mathematical  framework  for  quantifying  and  describing  the 
effects  of  these  low  dimensional  “organized  structures”  on 
Lagrangian  transport.  Consequently,  this  type  of  mathematical 
approach  should  prove  ideal  for  analyzing  transport  and  mix¬ 
ing  processes  associated  with  the  large  scale,  organized  mo¬ 
tions  observed  in  geophysical  flows. 

The  mathematical  development  of  dynamical  systems 
theory  is  also  being  influenced  by  studies  of  Lagrangian 
transport  problems  in  geophysical  fluid  dynamics.  The  flow 
studied  in  detail  in  this  article  was  two-dimensional  and  time- 
periodic.  Clearly  dynamical  systems  techniques  for  three-di¬ 
mensional  time-dependent  flows  must  be  developed  if  these 
methods  are  to  have  large  impact  on  Lagrangian  transport 
issues.  In  this  setting,  a  variety  of  new  mathematical  issues 
must  be  addressed.  The  beginnings  of  the  development  of  this 
theory  can  be  found  in  Beigie  et  al.  [1991],  [1992],  [1994], 
Mezic'  and  Wiggins  [1994],  and  Wiggins  [1992]. 

Finally,  we  mention  a  possibly  revolutionary  new  devel¬ 
opment  for  Lagrangian  transport  problems.  When  the  dynami¬ 
cal  systems  techniques  are  developed  so  that  they  can  be 
coupled  with  modem  computational  methods  for  computing 
velocity  fields  the  applicability  of  the  methods  will  be  greatly 
extended,  as  well  as  be  more  accessible  to  geophysical  fluid 
dynamicists.  Over  the  past  15  years  computational  fluid  dy¬ 
namics  has  developed  into  a  subject  in  its  own  right.  Now  we 
have  accurate  algorithms  for  solving  the  Navier  Stokes  equa¬ 
tions  in  a  variety  of  physically  important  settings.  However, 
many  problems  related  to  mixing  and  transport  begin  once  this 
step  has  been  accomplished.  That  is,  first  a  solution  to  the 
Navier  Stokes  equations  must  be  obtained  in  order  to  study  the 
transport  and  mixing  properties  associated  with  that  velocity 
field.  In  the  vast  majority  of  situations  arising  in  applications, 
this  solution  can  only  be  obtained  numerically.  Thus,  we  only 
have  a  numerical  representation  of  the  velocity  field.  Never¬ 
theless,  many  dynamical  systems  results  do  not  care  about  the 
specific  analytical  form  of  the  dynamical  system  under  con¬ 
sideration.  Rather,  they  require  that  only  certain  geometrical 
features  be  present.  For  example,  the  existence  of  stable  and 
unstable  manifolds  of  some  invariant  set  requires  only  the 
existence  of  a  hyperbolic  invariant  set  and  the  existence  of 
Smale  horseshoe  type  chaos  requires  only  the  transverse  inter¬ 
section  of  the  stable  and  unstable  manifolds  of  a  hyperbolic 


periodic  orbit.  Currently  our  research  is  focused  on  developing 
the  theory  and  algorithms  needed  to  perform  dynamical  sys¬ 
tems  type  analysis  with  numerical  representations  of  vector 
fields. 
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Editor’s  Note  -  Figures  4  through  10  show  various  design  arrangements  for  the  perturbation  of  fluid  flow  with  the  cusped 
jet  velocity  profiles.  Unfortunately,  the  red  and  blue  distinctions  indicated  in  the  captions  for  these  figures  are  not  shown 
However,  the  line  design,  line  breakup,  and  line  strength  are  discernable.  On  the  caver,  figures  8,  9,  and  10  of  this  article 
are  shown  with  color  distinctions. 


1.  Introduction 

The  difficulties  involved  in  collecting  data  on  oceanic 
motion  are  evident  and  the  inherent  uncontrollability  of  the 
ocean  for  experimentation  presents  a  formidable  scientific 
challenge  to  those  seeking  to  understand  it.  The  data  that  is 
collected  from  the  tracking  of  strategically  placed  floats  fits 
very  well  with  the  Lagrangian  perspective  on  fluid  mechan¬ 
ics,  wherein  motion  is  viewed  through  the  following  of  fluid 
particles,  as  it  is  assumed  that  the  paths  of  these  floats 
reflect  those  of  fluid  particles.  These  Lagrangian  trajecto¬ 
ries  are,  from  the  dynamical  systems  viewpoint,  the  orbits 
generated  by  a  vector  field  that  is  precisely  the  Eulerian 
velocity  field  of  the  fluid  flow.  This  perspective  thus  offers 
a  natural  mathematical  context  through  which  to  study 
questions  about  oceanographic  flows. 

However,  this  form  of  oceanographic  data  collection 
through  the  tracking  of  floats  suggests  a  deeper  reason  for 
the  suitability  of  the  dynamical  systems  perspective.  Since 
only  a  small  number  of  SOFAR  (or  RAFOS)  floats  can  be 
placed  in  areas  of  the  ocean  under  study,  and  they  can  only 
be  tracked  for  a  certain  period  of  time,  usually  a  year  or  less, 
the  resulting  particle  paths  must  be  viewed  as  a  sampling  of 
the  corresponding  Lagrangian  trajectories.  The  pictures  of 


these  float  paths  are  put  together  to  form  what  are  commonly 
known  as  “spaghetti” plots,  see  Figure  1  for  an  example  taken 
from  reference  4.  It  is  clear  that  we  can  only  rely  on  general 
qualitative  features  of  such  plots  for  an  assessment  of  the  fluid 
motion.  One  cannot  hope  to  give  an  explanation  of  detailed 
properties  of  these  paths.  Indeed,  such  an  explanation  would 
not  be  meaningful  as  any  individual  path  is  the  inevitable 
victim  of  small-scale,  temporary  effects.  However,  the  large- 
scale  features  of  the  paths  that  are  present  in  a  number  of 
samplings  should beexplained. 

The  modem  theory  of  dynamical  systems,  with  its 
emphasis  on  qualitative  properties  and  characteristic  geo¬ 
metric  features  of  trajectories,  provides  a  fitting  perspective 
on  the  issues  raised  by  the  study  of  oceanographic  flows. 
The  underlying  tenet  of  geometric  dynamical  systems  the¬ 
ory  is  that  certain  basic  geometric  structures  can  be  isolated 
which  organize  the  overall  structure  of  the  flow.  In  particu¬ 
lar,  certain  specific  trajectories,  if  they  can  be  located,  will 
demarcate  different  regions  of  qualitatively  distinctive  tra¬ 
jectory  motion.  Specifically,  stagnation  points  of  the  flow 
and  their  attendant  stable  and  unstable  manifolds  (sets  of 
trajectories  approaching  the  stagnation  point  in,  respec¬ 
tively,  forward  and  backward  time)  play  definitive  roles  in 
organizing  the  flow. 
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From  the  dynamicist’s  perspective,  the  oceanographic 
flow  fields  are  presented  to  us  on  three  levels  of  increasing 
accuracy,  and  complexity.  The  simplest  level  of  model  is  of  a 
velocity  field  that  is  postulated  as  an  accurate  model  of  the 
phenomenon  under  investigation.  This  will  be  formed,  com¬ 
monly,  in  terms  of  an  unstable,  or  neutral,  mode  superimposed 
on  a  base  flow.  Such  a  flow  leads  to  an  analytic  expression  for 
an  ordinary  differential  equation  whose  analysis  is,  in  princi¬ 
ple,  possible  by  existing  techniques  and  only  restricted  by  its 
complexity,  in  particular  its  dimension.  The  next  level  of 
model  is  the  numerical  model  which  consists  of  a  velocity  field 
given  as  a  numerical  database  that  is  constructed  by  solving 
the  relevant  partial  differential  equation  on  an  appropriate 
domain.  A  dynamical  systems  approach  to  such  a  velocity  field 
is  severely  impeded  by  the  fact  that  all  the  developed  tech¬ 
niques  are  based  on  the  idea  that  the  velocity  field  is  supplied 
to  the  mathematician  analytically.  A  new  perspective  is  thus 
called  for  and  the  development  of  techniques  for  handling 
numerically  generated  velocity  fields  is  a  significant  task  that 
remains  to  be  adequately  addressed.  In  particular,  the  numeri¬ 
cal  determination  of  the  relevant  structures,  such  as  stable  and 
unstable  manifolds,  for  numerical  velocity  field  data  sets  is 
required.  A  strong  case  has  been  made  by  Ridderinkhof  and 
Zimmerman^^,  for  the  relevance  of  these  dynamical  structures 
in  a  numerical  model  of  the  Wadden  Sea.  The  final  level  is  the 
experimental  data  set  available  for  the  situation  under  inves¬ 
tigation.  These  data  sets,  as  discussed  above,  are  usually  based 
on  the  tracking  of  floats  and,  while  insufficient  to  be  viewed 
as  a  complete  model,  must  always  be  used  as  the  final  arbiter 
of  the  success  of  a  given  piece  of  analysis. 

In  this  review,  we  shall  consider  models  at  the  first  level. 
The  formulation  of  such  models  as  base  flow  plus  perturbing 
mode  can  offer  great  insight  into  charcateristic  flows  that  may 
truly  exist  as  the  nonlinear  equilibration  of  such  “neutral”  or 


Figure  1. 

Spaghetti  diagram  ofRAFOS  floats  in  the  northern  Atlantic  Ocean 


SPAGHETTI  DIAGRAM  OF  RAFOS  FLOATS 
1984-1985 


“unstable”  directions  from  the  base  flow.  Explicit  analytic 
expressions,  or  even  approximations,  for  such  fully  nonlinear 
Eulerian  flow  fields  are  comparatively  rare,  and  these  “linear” 
approximations  offer,  with  the  existing  technology,  the  best 
way  to  make  meaningful,  analytic  predictions  about  the  nature 
of  Lagrangian  trajectory  motion  in  the  full  flow  field.  How¬ 
ever,  a  goal  of  the  mathematical  development  motivated  by 
oceanographic  problems  must  be  the  construction  of  better 
analytic  models  that  incorporate  real  nonlinear  effects  while 
remaining  faithful  to  the  underlying  physics. 

Our  focus  will  be  on  heteroclinic  orbits,  which  are  trajec¬ 
tories  joining  different  stagnation  points,  and  are  often  known 
as  separatrices  on  account  of  their  function  in  separating 
different  flow  behaviors.  We  favor  here  the  more  modem  term 
of  “heteroclinic  orbit”  because  the  main  point  of  the  message 
is  that,  under  perturbation,  these  structures  no  longer  perform 
this  function  of  separation.  Nevertheless,  there  remain,  typi¬ 
cally,  heteroclinic  orbits  after  perturbation.  We  shall  present 
three  examples  that  illustrate  how  an  elucidation  of  the  struc¬ 
ture  of  these  orbits,  in  particular  for  small  perturbations  of  the 
base  flow,  can  reveal  the  mechanisms  in  the  flow  for  the 
transport  of  particles. 

The  first  example  concerns  the  settling  of  a  particle  in  a 
cellular  flow  field.  The  particle  of  interest  here  is  not,  in  fact, 
a  fluid  particle  but  an  aerosol  particle,  and  a  key  point  in  the 
analysis  is  the  small  inertial  effect  experienced  by  such  a 
particle,  the  question  of  interest  is  whether  the  particle  under 
consideration  can  be  trapped  in  the  cell  or  continue  to  settle 
through  the  flow  indefinitely. 

An  interesting  feature  of  jets  in  the  ocean,  such  as  the  Gulf 
Stream,  is  the  spawning  of  rings  that  form  a  vortex  motion, 
and  this  is  the  motivation  for  our  second  example. 

These  rings  carry,  for  instance,  cold  water  into  the  region 
of  warm  water  to  the  south  of  the  Gulf  Stream.  They  tend  to 
travel  in  the  reverse  direction  from  the  jet  and  survive  for  a 
significant  amount  of  time,  often  until  hitting  the  continental 
shelf.  The  formation  of  these  rings  is  not  well  understood,  see 
Pratt  et  al.^®  and  many  of  their  features  remain  to  be  explained. 
Kirwan  et  al.^^’^^  formulated  a  model  for  such  rings  in  terms 
of  rotating  modons.  These  are  special  solutions  of  the  qua- 
sigeostrophic  equations  and  their  mixing  properties  will  be 
discussed  in  section  3. 

The  third  example  concerns  the  transport  of  water  within, 
and  across,  an  oceanic  jet.  The  general  mechanisms  for  mixing 
across  a  jet,  such  as  the  Gulf  Stream,  are  of  fundamental 
importance.  Various  kinematic  models  have  been  formulated 
for  the  clarification  of  this  transport,  but,  recently,  models  of 
velocity  fields  that  are  dynamically  consistent  have  been  de¬ 
rived  by  Pratt  et  al.^^  We  discuss  in  section  4  some  of  the  issues 
involved  in  the  transport  of  fluid  across  such  a  jet. 

In  each  of  these  examples,  there  is  an  “unperturbed”  flow 
field  thatis  aone-degree-of-freedom  Hamiltonian  system.  The 
distinguished  orbits  that  connect  stagnation  points  to  each 
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Figur*  2. 

A  diagrammatic  representation  of  the  streamlines  of  Stommel's 
model.  The  protected  pods  can  be  seen  in  each  of  the  four  cells 
depicted.  Trajectories  in  these  pods  correspond  to  particles  that 
are  forever  suspended,  while  the  particles  outside  will  settle. 


other  separate  regions  of  contrasting  Lagrangian  motion  from 
each  other.  The  perturbations  supplied  by  effects  unaccounted 
for  in  the  base  flow  will  tend  to  break  these  separating  orbits 
and  offer  the  possibility  of  transport  between  regions  that 
were,  in  the  unperturbed  case,  characterized  by  different  mo¬ 
tions.  We  shall  see,  in  each  of  the  examples  discussed  below, 
the  implications  of  the  perturbations  on  these  distinguished 
orbits. 

2.  Settling  in  a  Cellular  Flow  Field 

The  problem  of  whether  particles  of  non-fluid  matter 
settle  in  the  ocean  under  the  influence  of  gravity  is  obviously 
of  great  importance.  A  cellular  flow  field  is  a  model  of  Lang¬ 
muir  circulation  and  we  consider  here  the  issue  as  to  whether 
aerosol  particles  will  tend  to  settle  through  such  a  flow  field 
or  be  trapped  and  recirculate  within  that  cell.  A  key  point  is  the 
role  of  the  inertia  of  the  particle  in  question.  The  inertial  effects 
can  provide,  as  will  be  seen  below,  an  effective  dissipation  of 
the  Hamiltonian  that  arises  as  the  streamfunction  of  a  2-dimen¬ 
sional  incompressible  flow. 

The  settling  of  an  aerosol  particle  through  a  cellular  flow 
field  was  studied  by  StommeP^.  Stommel’s  analysis  neglected 
the  inertia  of  the  particle,  and  thus  the  particle  was  effectively 
treated  as  a  fluid  particle  settling  through  the  flow  field  under 
the  combined  influences  of  gravity  and  the  surrounding  fluid. 


He  showed  that,  in  this  scenario,  both  trapping  and  settling 
were  possible.  Indeed,  the  ratio  of  trapped  particles  to  settling 
particles  varied  with  the  characteristic  terminal  settling  veloc¬ 
ity  due  to  gravity.  In  physical  space,  when  the  influence  of 
inertia  is  omitted,  heteroclinic  orbits  to  the  stagnation  points 
on  vertical  invariant  lines  between  cells  separate  the  regions 
from  which  particles  will  settle  from  those  in  which  particles 
will  be  trapped,  as  shown  in  Figure  2. 

Indeed,  these  heteroclinics  form  invariant,  protected 
pods.  Under  the  periodic  velocity  field,  the  Lagrangian  trajec¬ 
tories  outside  a  pod  will  settle  through  all  of  the  cells  iey 
enter,  whereas  the  Lagrangian  trajectories  inside  a  pod  will  be 
trapped  forever. 

Despite  the  apparent  possibility  of  many  particles  being 
suspended  in  the  cells,  Maxey  and  Corrsin^^’^'*  observed,  based 
on  numerical  experiments,  that  inertial  particles  had  a  greater 
tendency  than  fluid  particles  to  settle.  The  true  problem  with 
inertia  included  is,  however,  a  perturbation  of  that  considered 
by  Stommel,  and  we  must  consider  whether  any  analogue  of 
the  protected  pod,  which  might  guaranteee  the  suspension  of 
particles,  exists  in  the  inertial  problem. 

If  the  inertia  of  the  particles  moving  through  the  cellular 
flow  field  is  small,  but  nonzero,  the  velocity  field  for  these 
particles  becomes  a  singularly  perturbed  version  of  that  for 
fluid  particles  acted  upon  by  the  cellular  flow  field  and  gravity. 
The  analysis  and  interpretation  of  such  singularly  perturbed 
problems  contains  many  pitfalls.  Fortunately,  a  technology 
exists  which  can  render  a  complete  geometric  description  of 
the  trajectories  of  such  perturted  systems,  and  it  can  be  in¬ 
voked  here  to  give  a  clear  and  unambiguous  answer  as  to 
whether,  under  the  influence  of  inertia,  particles  will  prefer  to 
settle  or  not. 

This  set  of  techniques  is  based  on  a  striking  collection  of 
theorems  due  to  Fenichel.^’^’^^  These  results  give  conditions 
under  which,  for  the  perturbed  system,  there  exists  an  invariant 
manifold  (a  smooth  surface  of  solutions)  that  carries  a  vector 
field  which  is  a  small  (regular)  perturbation  of  the  formal  limit 
equation. 

If  we  let  y(t)  denote  the  position  of  the  center  of  the 
particle  in  space  (y=(yj,y2))  and  take  the  downward  direction 
to  be  the  direction  of  positive  y2  (the  opposite  of  the  conven¬ 
tion  in  reference  1 1)  then  the  equations  of  the  following  model 
govern  the  motion  of  an  aerosol  particle  in  a  cellular  flow  field 

y  =  vi 

y  =  V2 

evj  =  _  vj  -I-  sin  cos  y2 
ev2  +  -  V2  -  w  -  cos  yj  sin  ^2  (1) 

Here,  v(t)  is  the  velocity  of  the  particle  and  W  is  the 
settling  velocity  scaled  so  that  0  <  W<  1.  The  small  parameter 
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Here,  v(t)  is  the  velocity  of  the  particle  and  W  is  the 
settling  velocity  scaled  so  that  0  <  W<  1.  The  small  parameter 
e  >  0  is  the  Stokes  number  and  measures  the  inertial  response 
time  of  the  particle  to  the  medium.  The  case  e  =  0  renders  the 
formal  limit  in  which  the  particles  have  no  inertia  and  behave 
merely  as  fluid  particles.  Fenichel’s  work  implies  that  there 
exists,  for  system  (1)  with  small  inertia,  an  invariant  manifold 
of  the  type  discussed  above.  This  invariant  manifold,  known 
as  the  slow  manifold,  is  given  by  a  small  perturbation  of  the 
algebraic  conditions  found  by  formally  setting  e  =  0  in  the  last 
two  equations  of  (1).  An  application  of  Fenichel’s  further 
results  combined  with  some  simple  dynamical  systems  argu¬ 
ments  show  that,  in  this  case,  the  slow  manifold  is  actually 
globally  attracting.  Thus  the  asymptotic  motion  of  all  La- 
grangian  trajectories  will  be  determined  by  the  vector  field  on 
the  slow  manifold. 

A  key  point  to  notice  is  that  the  slow  manifold  is  para¬ 
metrized  by  the  physical  space  variables  (yi,y2)-  To  analyze 
the  true  behavior  on  the  manifold,  it  is  necessary  to  calculate 
the  0(e)  terms  in  the  perturbed  vector  field.  This  is  achieved 
by  using  the  invariance  property  of  the  manifold  and  the  full 
equation  itself.  The  resulting  system  on  the  slow  manifold  is 

jij  =  sin  yj  cos  >2  ~  ^  sin  >1  0^  si"  ^i)  +  ^ 

(2) 

y^=i  —  W-  COS sin y2  +  E  cos  y2  (W cos +  sin y2)  +  )• 

As  commented  above,  the  long-time  behavior  of  aerosol 
particle  trajectories  is  determined  completely  by  the  trajecto¬ 
ries  of  this  vector  field.  Moreover,  the  interpretation  of  these 
trajectories  is  considerably  facilitated  by  the  fact  that  the 
dependent  variables  in  (2)  are  exactly  the  space  variables. 
Setting  e  =  0  in  (2),  we  recover  the  velocity  field  for  fluid 
particles,  the  trajectories  of  which  are  sketched  in  Figure  2. 
We  are  now  ready  to  address  the  key  question  as  to  whether 
particles  should  be  expected  to  settle  under  the  influence  of 
the  vector  field  (2).  The  answer  will  obviously  come  from 
studying  the  fate  of  the  pods  present  in  the  e  =  0  system  under 
the  perturbation  supplied  by  the  0(e)  terms  of  (2).  This  is 
tantamount  to  asking  how  the  heteroclinic  orbit  surrounding 
the  pod  in  the  interior  of  each  cell  will  break  under  the 
influence  of  this  perturbation  (the  vertical  heteroclinic  stays  in 
place).  A  few  simple  sketches  will  convince  the  reader  that  if 
the  unstable  manifold  falls  outside  the  stable  manifold  then  the 
process  of  settling  is  facilitated,  whereas  the  reverse  configu¬ 
ration  would  result  in  more  trapping.  The  determination  of  the 
direction  of  breaking  of  heteroclinic  orbits  is  the  objective  of 
Melnikov  calculations,*®’^®  and  the  issue  can  be  resolved  by 
calculating  the  sign  of  a  certain  integral.  This  integral  presents 
many  technical  difficulties,  and  we  do  not  evaluate  it  directly, 
but  it  can  be  shown  analytically  that  its  actual  sign  corresponds 


to  the  unstable  manifold  falling  outside  the  stable  manifold. 
This  indicates  a  strong  enhancement  of  settling". 

In  this  way,  we  show  that  any  small  inertial  effects  in  fact 
force  particles  to  setUe  through  a  cellular  flow  field  and  make 
the  probability  of  suspension  almost  non-existent. 

Figure  3  provides  a  sketch  of  the  resulting  flow  of  parti¬ 
cles  with  small  nonzero  inertia. 

It  is  also  shown  analytically  in  reference  11  that  certain 
characteristic  paths  emerge  in  the  slow  manifold,  along  which 
the  aerosol  particles  will,  asymptotically,  congregate  and 
travel.  This  corroborates  further  numerical  observations  due 
to  Maxey  and  Corrsin.*^’*'* 

Although  the  underlying  cellular  flow  field  is  a  reasonable 
model  in  certain  oceanographic  circumstances,  it  is,  of  course, 
very  idealized.  Nevertheless,  our  work  shows  that  the  inclu¬ 
sion  of  inertial  effects,  if  properly  interpreted  on  a  slow  mani¬ 
fold,  supplies  an  effective  dissipation  of  the  Hamiltonian 
streamfunction.  Such  effects  are  also  seen  in  recent  related 
work  of  Tio  et  al.^*,  in  which  settling  through  a  critical  layer 
is  studied.This  may  partly  account  for  observations  in  which 
particulate  matter  in  the  ocean  is  seen  to  collect  in  certain 
well-defined  areas^®’^^,  despite  the  fact  that  it  is  under  the 
influence  of  a  time-dependent,  incompressible,  almost  strati¬ 
fied  flow,  which  one  would  otherwise  expect  to  be  srtongly 
homogenizing. 

3.  Rotating  Modons 

In  recent  work,  Mied,  Kirwan  and  Lindemann*^  have 
proposed  rotating  modons  as  models  of  rings  that  result  from 


Figure  3. 

The  trajectories  of  those  particles  lying  on  the  slow  manifold  are 
depictedhere.  Other  trajectories  will  asymptotically  approach  this 
configuration.  Small  inertia  is  included  and  has  the  effect  of 
breaking  the  heteroclinic  orbits  in  such  a  way  as  to  force  settling. 
Particle  paths  do  not  become  trapped  and  those  lying  within  the 
protected  pods  in  the  zero  inertia  case  can  leak  out  and  settle. 
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Figure  4. 

Hamiltonian  (effective  streamfunction)  of  the  feature  model  for 
82-B.  Solid  and  dashed  curves  are  positive  and  negative  values, 
respectively.  The  contours  have  been  normalized  by  a  convenient 
plotting  scale.  The  dotted  circles  show  model  scale  lengths  for 
the  upper  and  lower  model  layers.  The  inner  circle  radius  is  120 
km.  There  are  two  types  of  critical  points  in  this  plot  The  three 
diamonds  depict  stagnation  points  that  are  local  extrema.  There 
are  two  dots  near  the  tips  of  the  SP  contour  lines.  These  are  saddle 
points.  The  SP  contour  is  the  value  of  teh  Hamiltonian  at  the 
saddle  points. 


the  pinching  off  of  vortex-like  regions  from  such  jets  as  the 
Gulf  Stream.  These  multipole  vortex  systems  consisting  of  one 
or  more  cyclonic  eddies  moving  around  an  anticycionic  ring 
often  are  seen  in  AVHRR  imagery  in  the  coastal/open  ocean 
interface.  The  nonlinear  feature  model  developed  by  Kirwan 
et  al.^^  accounts  for  many  of  the  charcteristics  seen  in  the 
imagery.  These  modons  are  specific  solutions  of  the  qua- 
sigeostrophic  potential  vorticity  equations  for  a  2-layer  model, 
in  such  a  setting  they  were  first  introduced  by  Flierl  et  al.^. 
Modons  are  steady  in  a  frame  that  is  rotating  at  a  fixed 
frequency.  Bottom  topography  can  be  included,  and  this  can 


have  the  effect  of  trapping  the  modon  and  preventing  it  from 
lateral  movement,  see^*.  Of  particular  interest  are  dipolar 
vortices  in  which  two  regions  of  opposite  rotation  are  matched 
together.  It  is  suggested  in  reference  16  that  such  dipoles  model 
pairs  of  warm  and  cold  cores  that  are  attached  to  each  other, 
and  can  even  become  so  attached  as  the  result  of  a  collision. 
Of  interest  in  our  study  is  the  local  mixing  properties  of  such 
modons.  For  instance,  is  it  possible  that  fluid  from  outside  the 
modon  can  be  mixed  into  it?  Can  it  flow  into  the  core?  Or 
between  the  two  complementary  cores?  If  not,  what  external 
mechanisms  might  facilitate  such  mixing? 

The  rotating  modons  are  solutions  of  the  quasigeostrophic 
equations  of  a  prescribed  form.  They  are  constructed  by  find¬ 
ing  solutions  in  different  circular  regions  and  matching  at  the 
boundaries.  This  is  quite  a  complicated  procedure  that  leads  to 
formidable  calculations  for  solutions  of  the  resulting  partial 
differential  equations. 

In  the  paper  by  Kirwan  et  al.  each  layer  is  divided  into 
an  inner  (circular)  modon  region  and  an  outer  region.  The  radii 
of  the  regions  in  each  layer  will,  in  general,  be  different.  The 
quasigeostrophic  equations  for  the  potential  vorticity  q.  and 
the  streamfunction  \|/.  are  assumed  to  hold  in  the  regions  of 
each  layer,  i  =  1 ,2,  with  appropriate  matching  conditions  at  the 
boundary.  Solutions  of  these  equations  which  are  steady  in  a 
rotating  frame  of  frequency  co  can  be  found  under  the  ansatz 
that  the  effective  streamfunction 


Figure  5. 

This  vortex  system  has  a  rotation  period  of  30.6  inertial  days. 
Positive  values  are  depicted  as  solid  curves  and  negative 
values  as  dashed  curves.  Two  blobs  are  initialized  on  either 
side  of  a  saddle  point;  the  red  blob  is  just  inside  the  saddle 
point  while  the  blue  blob  is  just  outside.  The  blob  radii  are  5km 
and  there  is  a  2  km  distance  between  the  blobs. 
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Figure  6. 

Flow  associated  to  the  single  cusp  jet  model  with  a  single  neutral 
mode  superimposed.  When  viewed  in  a  reference  frame  moving 
with  the  speed  of  the  neutral  wave,  the  flow  is  time-independent 
and  admits  two  rows  of  critical  points  at  the  steering  lines.  Fluid 
between  the  steering  lines  moves  to  the  right  as  its  speed  is 
greater  than  the  speed  of  the  reference  frame,  and  fluid 
outside  the  steering  lines  moves  to  the  left.  In  each  of  the 
figures  of  this  section,  the  heteroclinic  orbits  are  colored  red 
(blue)  for  the  pieces  that  are  generated  as  unstable  (stable) 
manifolds.  This  is  not  so  meaningful  in  the  unperturbed  case 
(Figures  6,  7  and  9)  since,  on  the  heteroclinic  orbits  the  stable 
and  unstable  manifolds  coinicide.  the  reamining  trajectories 
are  colored  red. 


X 


and  the  potential  vorticity  are  functionaliy  related,  with  a 
different  relation  in  each  layer  and  region.  The  quasigeostro- 
phic  equations  are  then 

J{Yi,q-)  =  Si, 

where  Sj  are  sources  at  the  boundaries  of  the  regions  to 
compensate  for  the  jumps  in  potential  vorticity.  The  relevant 
solutions  are  then  found  by  Kirwan  et  al.'*  from  solving  the 
partial  differential  equations  coming  from  the  functional  rela¬ 
tions. 

In  physical  variables,  namely  r,  0  and  t,  the  Lagrangian 
trajectories  are  those  of  a  Hamiltonian  system  with  time-de¬ 
pendent  forcing.  There  are  two  such  systems,  one  for  each 


layer,  and  the  Lagrangian  trajectories  in  each  layer  are  inde¬ 
pendent  of  those  in  the  other  layer,  although  the  corresponding 
streamfunctions  are  not.  Within  each  layer  the  Lagrangian 
trajectories  are  the  same  for  any  depth.  The  trajectory  equa¬ 
tions  are  then 

dr  _  1 

dt  r  ae 

(3) 

dt  r  dr  ' 

where  the  streamfunction  Vi  ^  '  ^0-  The  streamfunction 

will  be  continuous  across  the  boundaries  of  the  various  re* 
gions,  but  not  necessarily  smooth.  Since  the  system  (3)  is  a 
time  dependent  Hamiltonian  system,  its  trajectories  can  be, 
ostensibly,  chaotic.  However,  if  viewed  in  a  rotating  frame,  it 
is  easily  seen  that  the  effective  Hamiltonian  becomes  a  time 
independent  Hamiltonian  and  after  the  appropriate  changes  of 
coordinates,  (3)  can  be  recast  as  a  one-degree-of- freedom 
Hamiltonian  system.  We  introduce  the  variables  0  =  0  -  cot  and 
a  scaled  time  t,  so  that  (3)  becomes 


Figure  7. 

Increasing  the  wavenumber  of  the  neutral  mode  draws  the 
steering  lines  closer  to  the  centerline  of  the  jet  and  the  meander 
of  the  jet  becomes  more  pronounced.  The  assymetry  in  the 
cat's  eyes  is  due  to  the  exponential  form  of  the  eigenfunctions. 
Here  the  neutral  mode  has  wavenumber  k  =  8  and  amplitude 
a  =  0.01. 


X 
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Figure  8. 


Adding  a  second  neutral  mode  to  the  jet  model  results  in  a  time 
periodic  vector  field.  Here  the  second  mode  is  moving  faster 
than  the  first  mode  with  wavenumber  k2  =  16  and  amplitude  e 
=  0.001.  In  the  presence  of  this  time  periodic  perturbation  the 
heteroclinic  orbits  break  up  into  transverse  intersections  of  the 
stable  and  unstable  manifolds  for  the  associated  Poincare 
map.  The  most  significant  mixing  occurs  on  the  inside  of  the 
jet,  however,  no  mixing  occurs  across  the  body  of  the  jet.  The 
color  coding  (red-unstable,  blue-stable  manifolds)  shows 
clearly  the  intermingling  between  stable  and  unstable  mani¬ 
folds  on  the  inner  heteroclinic  orbits. 


dr 

di~  d@ 

(4) 

dl  dr 

where  the  right  hand  side  depends  only  on  the  phase  space 
variables  r  and  ©. 

This  casting  of  the  problem  has  rather  striking  implica¬ 
tions.  If  the  modon  is  viewed  in  the  rotating  frame,  then  the 
trajectories  must  follow  the  contours  of  the  effective  Hamil¬ 
tonian  and  are,  in  particular,  quite  tame.  We  shall  focus  on  the 
case  of  a  dipolar  modon.  In  recent,  and  ongoing,  work  Kirwan 
and  the  authors  have  discovered  stagnation  points  in  the  flow 
field  (3)  for  the  upper  layer.  These  correspond  to  periodic 


Lagrangian  trajectories  in  physical  space.  These  stagnation 
points  are  crucial  in  understanding  the  dynamics  of  (3).  Indeed 
they  are  located  near  the  edge  of  the  modon  boundary  (radius 
of  inner  region)  for  the  upper  layer.  The  outer  ring  of  hetero¬ 
clinic  orbits  (called  a  heteroclinic  cycle)  surrounds  both  the 
cyclonic  and  anti-cyclonic  vortices  (warm  and  cold  cores)  of 
the  modon.  The  inner  cycle  surrounds  only  the  cyclonic  vor¬ 
tex.  Figure  4  shows  the  qualitative  nature  of  the  phase  portrait 
of  (3).  The  figure  is  based  on  data  for  a  model  of  82-B,  using 
hydrographic  and  AVHRR  data.  This  data  is  well  documented 
in  Kenelly  et  al.^^  The  two  stagnation  points  are  located  in  the 
north-east  and  north-west  sides  of  the  modon  structure.  The 
heteroclinic  orbits  are  contours  of  the  effective  Hamiltonian 
joining  the  stagnation  points.  These  are  not  cpatured  exactly 
in  the  Figure  4,  but  reasonable  approximations  can  be  seen  in 
the  contours  marked  SP,  it  is  left  to  the  readers  imagination  to 
see  how  these  are  connected  to  the  saddle  stagnation  points  to 
form  the  heteroclinic  orbits.  Of  particular  significance  in  this 
analysis  is  that  the  heteroclinic  cycles  provide  a  barrier  to  the 
mixing  of  fluid  in  certain  regions.  The  outer  homoclinic  cycle 
prevents  any  outside  fluid  from  infiltrating  the  modon  struc¬ 
ture.  This  is  not  obvious  from  pictures  in  the  physical  space  as 
this  homoclinic  orbit  will  not  be  circular  and  hence  will  change 
position  as  the  modon  rotates.  Nevertheless,  if  a  fluid  particle 
is  outside  the  rotated  version  of  this  cycle  at  an  appropriate 
time,  the  modon  flow  will  never  carry  it  through  the  modon. 
In  some  sense,  this  cycle  defines  the  modon  as  it  supplies  the 
modon  boundary.  The  inner  heteroclinic  cycle  provides  the 
barrier  between  the  two  main  vortices  of  the  modon.  It  should 
be  noted  that  the  contouring  of  the  effective  streamfunction 
(Hamiltonian)  has  uncovered  a  third  vortex  patch  which  was 
somewhat  unexpected.  It  is,  however,  viewed  as  less  signifi¬ 
cant  than  the  two  primary  vortices. 

A  numerical  experiment,  performed  by  B.Lipphart  (Old 
Dominion  University),  gives  clear  corroboration  of  the  protec¬ 
tion  of  the  modon  by  the  outer  heteroclinic  cycle.  In  this 
numerical  experiment,  shown  in  Figure  5,  two  blobs  of  initial 
data  are  placed  in  locations  near  the  northeasterly  stagnation 
point.  The  outer  blob  (blue)  is  placed  just  outside  the  outer 
heteroclinic  cycle  at  the  stagnation  point,  and  the  inner  blob 
(red)  just  inside.  The  resulting  particle  paths  are  then  strobed 
at  six,  equally  spaced,  time  intervals  in  physical  space  ending 
after  one  period.  The  contour  lines  sketched  in  Figure  5  are 
those  of  the  original  streamfunction  (\|/^  and  not  Yj)  and  their 
presence  is  used  to  facilitate  the  viewing  of  the  modon.  The 
different  fates  of  the  two  blobs  can  be  seen  quite  clearly  from 
the  numerical  plots.  Indeed,  the  blue  particle  paths  are  kept  out 
of  the  modon  structure,  while  the  red  ones  are  sucked  into  the 
region  between  the  vortices.  This  difference  is  especially 
striking  when  one  considers  that  the  initial  separation  of  the 
blobs  is  not  that  great,  only  about  2km.  in  physical  space. 

In  the  real  oceanographic  systems,  one  cannot  expect  the 
modons  to  be  totally  protected  from  the  outside  fluid.  Indeed, 
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the  outside  fluid  will  be  sucked  through  the  center  of  the 
modon  and,  to  some  extent,  into  the  edges  of  the  vortices.  The 
question  then  arises  as  to  what  perturbing  mechanisms  might 
effect  such  mixing  of  fluid  through  the  modon  structure.  This 
is  the  subject  of  our  ongoing  work  in  collaboration  with 
Kirwan,  but  some  observations  can  be  made  as  to  why  and 
how  we  anticipate  this  being  realized.  The  outside  ocean  is  a 
ready  supplier  of  perturbations  that  will  manifest  themselves 
as  time  varying  influences  on  the  modon  structure.  It  is  a 
well-known  fact  in  dynamical  systems  that  for  area-preserving 
systems,  homoclinic  orbits  for  the  associated  Poincare  map 
will  persist  under  perturbation.  This  can  be  applied  here  but, 
since  the  perturbed  system  will  be  time-dependent,  has  the 
consequence  that  the  stable  and  unstable  manifolds  of  a  fixed 
point,  close  to  the  original  stagnation  point,  will  intersect. 
Since  these  are  now  invariant  manifolds  for  a  map,  their 
intersection  produces  a  quite  different  scenario.  Indeed,  the 
manifolds  for  the  perturbed  system  should  intersect  trans¬ 
versely,  which  we  expect  as  a  consequence  of  a  symmetry 
argument  similar  to  that  noted  in  the  following  section,  and 
chaotic  mixing  of  fluid  particles  between  the  regions  that  were 
previously  protected  from  each  other  is  now  likely,^*.  It  may 
be  commented  here  that  the  conservation  of  potential  vorticity 
will  prevent  chaotic  mixing,  as  observed  recently  by  Brown 
and  Samelson^.  However,  in  the  unperturbed  case,  the  poten- 

Figure  9. 

Modeling  the  dynamics  at  a  steering  line  where 
U'(y)  =  0.  Here  the  steering  line  occurs  at  a  point  where  the 
base  velocity  profile  attains  a  minimum  (W  (0)  0).  Because  of 
the  scaling  in  this  analysis  the  dynamics  actually  occur  over  a 
smaller  range  ofy  than  in  the  previous  figures. 


tial  vorticity  and  effective  streamfunction  are  linearly  depend¬ 
ent,  that  is  indeed  how  they  were  found  for  the  modon,  and 
thus  the  analysis  of  reference  5  does  not  apply  in  this  case. 
Even  though  we  expect  mixing  through  the  modon  to  occur 
under  such  time-dependent  perturbations,  the  KAM  Theorem, 
see  reference  1  for  instance,  can  be  applied  to  show  that,  under 
small  perturbations,  some  of  the  orbits  surrounding  the  vortex 
regions  will  survive  as  invariant  tori.  This  will  protect  the  vortices 
themselves  from  being  sucked  into  the  mixing  process.  Thus,  we 
expect  for  the  flow  that  fluid  will  be  mixed  through  the  modon, 
that  is  to  say  between  the  two  vortex  cores,  but  not  into  the  cores 
themselves. 

A  key  point  here  is  how  to  model  the  time  dependent 
perturbation,  whether  by  a  superimposed  mode  or  a  time 
variation  of  parameters  in  the  problem,  such  as  the  radii  of  the 
regions  in  one  of  the  layers.  Issues  of  dynamical  consistency 
and  physical  relevance  are  of  central  importance.  However, 
we  expect  that  mixing  of  fluid  will  occur  through  the  modon 
for  almost  all  such  perturbations.  Given  the  formulation  of  an 
appropriate  model,  mixing  rates  can  even  be  calculated  from 
the  Melnikov  function,  as  formulated  by  Rom-Kedar  and 
Wiggins^'^\  and  for  some  more  recent  theory^,  and  the 
different  mechanisms  available  can  be  compared  for  the  effi¬ 
cacy  in  the  process  of  mixing  through  the  modon. 

4.  Mixing  across  Jets 

The  mechanism  for  the  transport  of  fluid  particles  across 
such  oceanic  jets  as  the  Gulf  Stream  is  not  well  understood. 
Near  the  surface  the  presence  of  high  potential  vorticity  gra¬ 
dients  at  the  center  of  the  jet  appear  to  act  as  barriers  to  this 
transport.  Indeed,  this  has  been  experimentally  corroborated 
by  Behringer,  et  al.^.  Nevertheless,  it  appears  to  occur  either 
by  ring  detachment  or  straightforward  advection. 

Many  models  have  been  developed  with  respect  to  which 
these  questions  can  be  posed.  Bowei^  and  Samelson  studied 
kinematic  models  based  on  the  meander  structure  of  the  Gulf 
Stream  evident  from  spaghetti  plots,  such  as  shown  in  Figure 
1.  More  recently,  dynamically  consistent  models  were  formu¬ 
lated  by  del-Castillo-Negrete  and  Morrison^  and  Pratt,  Lozier 
and  Beliakova^'.  The  first  group  used  the  Bickley  jet  (sech^y 
profile)  as  the  base  flow  coupled  with  superimposed  neutral 
modes.  Pratt  et  al  used  the  cusped  jets  derived  from  a  1  Vi  layer 
model  in  the  paper  by  Pratt  and  Stem^^.  These  cusped  jets  have 
piecewise  constant  potential  vorticity,  and  the  center  of  the  jet 
is  characterized  by  a  discontinuity  in  the  potential  vorticity. 
Since  these  cusped  jets  are  linearly  stable,  it  is  reasonable  to 
model  perturbations  in  the  full  flow  field  by  superimposing  on 
the  base  flow  selected  neutral  modes.  It  is  assumed  that  two 
waves  are  used  and  the  first  is  a  larger  perturbation  than  the 
second.  The  first  neutral  mode  introduces  a  characteristic 
speed  and  the  steering  lines  occur  where  this  matches  the  speed 
of  the  fluid  particles  in  the  base  jet.  As  pointed  out  by  both 
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Figure  10. 

Adding  a  second  neutral  mode  with  a  speed  slightly  greater  than 
the  first  mode  results  in  significant  mixing  across  the  region  of  the 
steering  line.  This  is  evident  from  the  strong  breakup  of  the  stable 
(blue)  and  unstable  (red)  manifolds  of  the  associated  Poincare 
map.  Again  the  intermingling  is  shown  by  the  color  coding,  see 
Figure  7. 


groups,  in  the  regions  near  the  steering  lines  chaotic  trajecto¬ 
ries  can  appear  and,  more  importantly,  advective  mixing  be 
facilitated.  Such  chaotic  structures  will  be  expected  to  appear 
when  a  second  wave  is  superimposed  on  the  flow  field.  If  its 
wavenumber  and  frequency  are  unmatched  to  those  of  the 
primary  wave,  the  resulting  forcing  on  the  velocity  field  is  time 
dependent  and  such  a  scenario  is  known  to  produce  chaotic 
motion. 

A  comment  about  potential  vorticity  is  in  order  here,  see 
the  book  by  Pedlosky  for  the  relevant  definitions  and  discus¬ 
sion.  As  observed  recently  in  the  paper  by  Brown  and  Samel- 
son^  a  two-dimensional  incompressible  flow  with 
time-periodic  forcing,  which  preserves  potential  vorticity, 
cannot  be  chaotic  if  the  gradients  of  the  potential  vorticity  and 
the  streamfunction  are  independent.  Indeed,  the  Lagrangian 
trajectories  are  constrained  in  this  case  to  2-dimensional  sur¬ 
faces,  determined  by  the  potential  vorticity  contours,  and  will 
thus  not  be  chaotic.  As  commented  by  del-Castillo-Negrete 
and  Morrison^,  this  apparent  contradiction  is  resolved  by  the 
fact  that  flows  under  consideration  are  not  exactly  potential 
vorticity  conserving.  However,  from  the  mathematical  view¬ 
point,  this  breaking  of  the  potential  vorticity  conservation 
arises  from  the  approximations  made  in  implementing  the 


model.  These  approximations  are  made  in  assuming  that  the 
flow  field  is  given  as  a  base  flow  with  two  neutral  modes 
added.  The  full  nonlinear  field,  which  is  only  modelled  by  this 
linear  superposition,  is  potential  vorticity  conserving,  but  this 
may  not  hold  for  the  flow  field  that  is  actually  used.  Moreover, 
potential  vorticity  can  be  broken  by  the  numerical  approxima¬ 
tions  made  in  calculating  the  Lagrangian  trajectories.  Of 
course,  in  the  real  ocean  potential  vorticity  is  not  conserved 
and  we  do  expect  orbits  heteroclinic  to  stagnation  points  to 
break  up  in  such  a  way  as  to  produce  chaotic  motion.  It  would 
be  preferable,  however,  if  this  effect  was  modelled  directly 
rather  than  as  a  by-product  of  the  approximations  made.  In 
other  words,  we  would  hope  to  be  able  to  produce  these  effects 
of  chaotic  mixing  by  producing  models  that  more  accurately 
reflect  the  actual  ways  in  which  potential  vorticity  is  broken, 
rather  than  through  less  accurate  ways  of  modelling  systems 
in  which  it  is  conserved.  Nevertheless,  the  presence  of  chaotic 
mixing  is  credible  whenever  such  heteroclinic  orbits  appear 
and  we  believe  that  it  is  a  significant  mechanism. 

The  models  under  consideration  are  effectively  two-di¬ 
mensional,  and  the  effects  of  vertical  motion  may  be  quite 
significant.  It  has  been  suggested  that  mixing  across  a  jet  may 
occur  more  readily  at  deeper  levels,  see  Lozier  and  Riser^^. 
Combining  this  feature  with  the  possibility  of  vertical  motion 
offers  some  very  tantalizing  possibilities.  The  effects  of  depth 
can  be  studied  within  the  stratified  models  discussed  above  as 
the  depth  does  occur  as  a  parameter  within  the  flow  field. 
Questions  can  then  be  asked  about  the  dependence  of  advec¬ 
tive  mixing  on  depth,  for  instance:  is  mixing  facilitated  at 
lower  depths?  The  picture  one  should  keep  in  mind  is  that 
described  by  Pratt,  Lozier  and  Beliakova^^  in  which  the  base 
jet  is  bracketed  by  two  steering  lines  determined  by  the  pri¬ 
mary  wave.  On  the  surface,  these  steering  lines  are  well 
separated,  but  as  one  looks  deeper,  the  steering  lines  move 
toward  each  other  until  they  converge  and  cancel  each  other, 
at  a  certain  depth.  Below  that  depth  the  fluid  panicles  of  the 
base  jet  all  move  slower  than  the  primary  wave. 

To  study  the  dynamics  of  the  Lagrangian  trajectories  in 
the  vicinity  of  steering  lines  we  began  with  the  analytical 
model  for  a  cusped  jet  velocity  profile  formulated  by  Pratt, 
Lozier  and  Beliakova^^  The  flow  is  taken  to  be  two  dimen¬ 
sional  and  is  governed  by  a  streamfunction  \|/(x,y,t).  The 
Lagrangian  particle  motions  satisfy  the  equation 

dx  _d\y 

dt  dy 

dy  3\|/ 

di  dx 

The  base  flow  of  the  jet  is  a  shear  flow  with  velocity  profile 

U(y)  =  e-^y\ 
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=  aOCy)  cos  k(x  -  ct), 
where  the  eigenfunction  <I)(y)  is  given  by 
OCy)  = 

The  wavenumber  k,  the  wave  speed  c,  and  the  cross-stream 
decay  rate  k  are  related  by 

In  a  frame  moving  with  speed  c,  the  Lagrangian  trajecto¬ 
ries  are  solutions  of  the  equations 

^  =  [/ (y)  -  c  -  a  O'O’)  cos  A:C. 


^  =  -ak0(y)  sin^C 

at 

where  we  have  introduced  the  variable  ^  =  x  >  ct.  As  shown  in 
this  system  will  have  two  steering  lines,  for  appropriate 
values  of  c,  and  the  structure  of  the  trajectories  near  these  lines 
will  be  of  a  row  of  cat’s  eyes,  familiar  from  the  Stuart  vortex 
flow,  see  Figure  6.  Increasing  the  speed  c  of  the  superimposed 
wave  draws  the  steering  lines  closer  to  the  centerline  of  the  Jet. 
The  exponential  form  of  0(y)  introduces  asymmetry  into  the 
shape  of  the  periodic  orbits  and  the  meandering  shape  of  the 
jet  is  more  pronounced,  as  seen  in  Figure  7. 

The  break-up  of  these  heteroclinic  trajectories  is  forced 
by  a  second  neutral  wave  being  superimposed  on  the  stream- 
function.  This  secondary  wave  is  of  the  same  form  as  the  first 
wave, 

\|;2  (x:,y,0  =  £02(y)  cos  kj  (x  -  C2O. 

but  the  wave  speed  C2  is  chosen  different  from  the  speed  c  of 
the  first  mode.  The  wavenumber  k2  and  speed  C2  are  related  as 
in  (5),  and  the  amplitude  e  is  taken  to  be  small  relative  to  a. 
The  Lagrangian  trajectories  are  now  governed  by  a  time-de¬ 
pendent  vector  field  where  the  order  t  time  dependency  is 
periodic  with  angular  frequency  O)  =  k2(c2'C).  The  vector  field 
is  of  the  form 

^^U(y)-c-a<S>'  (y)  cos - eO'2 (y)  cos  (^2^ - “0- 


^  =  -alc^(y)  smkC,-£k02{y)  sin  {k^^-cat)- 
at 

Since  (7)  is  a  time-dependent  Hamiltonian  system  for  e 
0,  the  secondary  wave  is  expected  to  break  the  heteroclinic 


orbits  on  either  side  of  each  set  of  cat’s  eyes.  This  can  be 
verified  by  calculating  the  Melnikov  function  associated  with 
the  problem.  This  function  gives  the  leading  order  term  in  the 
signed  distance  between  the  stable  and  unstable  manifolds  of 
two  unstable  periodic  solutions  which  are  created  by  the 
perturbation  in  the  vicinities  of  the  two  unperturbed  stagnation 
points.  For  the  vector  field  (7)  the  Melnikov  function  can  be 
written  in  the  form 

where  x'’(t-to)  is  the  unperturbed  heteroclinic  orbit  connecting 
the  stagnation  points  for  e  =  0,  and  {.,.)  denotes  the  Poisson 
bracket  in  the  ^  and  y  variables.  The  explicit  form  of  x'’(t-tQ) 
seems  hopeless  to  determine,  but  it  is  in  fact  not  needed  for 
finding  zeros  of  the  Melnikov  function.  It  turns  out  that,  using 
the  symmetry  of  x'’(t-t(,)  with  respect  to  the  fixed  C  line  through 
the  center  of  the  unperturbed  problem,  one  can  argue  that  M(to) 
admits  isolated  zeros.  A  similar  reasoning  yields  that  these 
zeros  are  transverse,  therefore  the  heteroclinic  connection 
between  the  unperturbed  stagnation  points  breaks  up,  but 
individual  transverse  heteroclinic  solutions  between  nearby 
unstable  periodic  motions  survive. 

This  break-up  can  be  observed  numerically  through  the 
distinctive  looping  caused  by  transverse  intersections  and  the 
associated  heteroclinic  tangle.  Figure  8  shows  the  result  of 
perturbing  the  flow  shown  in  Figure  7  for  the  case  when  the 
wavenumber  k2  of  the  secondary  mode  is  twice  the  wavenum¬ 
ber  of  the  first  neutral  mode.  An  interesting  feature  here  is  that 
it  is  significantly  harder  to  detect  the  presence  of  transverse 
intersections  along  the  outside  edge  of  the  jet.  This  is  similar 
to  a  feature  noted  by  Samelson^  for  the  kinematic  meander 
model.  The  interesctions  on  the  outer  edge  seem  to  occur  at  a 
much  higher  order,  and  it  appears  that  the  outside  heteroclinics 
do  not  break  at  all.  Whereas,  when  the  steering  lines  draw 
closer  to  one  another,  the  heteroclinics  in  the  interior  of  the  jet 
are  easily  detected.  In  agreement  with  the  analysis  for  the 
derivative  of  the  Melnikov  function,  we  have  observed  nu¬ 
merically  that  the  strong  heteroclinic  breakup  can  be  made  to 
occur  on  the  opposite  side  of  the  loop  by  reversing  the  sign  of 
the  angular  frequency  (o.  This  implies  that  the  regions  of 
significant  mixing  may  depend  on  the  relative  speeds  of  the 
two  superimposed  neutral  waves,  since  the  sign  of  (o  depends 
on  the  difference  of  the  wave  speeds  c  and  C2.  In  particular,  if 
the  secondary  mode  is  faster  than  the  primary  neutral  mode, 
then  the  resulting  mixing  is  more  intense  within  the  body  of 
the  jet,  whereas  the  boundary  of  the  body  of  the  jet,  viewed  as 
being  the  two  steering  lines,  is  virtually  intact.  As  a  result, 
fluids  particles  are  kept  outside  of  the  jet  by  the  outer  stable 
and  unstable  manifolds  which  form  barriers  to  mixing.  At  the 
same  time,  if  the  wave  speed  of  the  secondary  mode  is  smaller 
than  that  of  the  primary  mode,  this  boundary  breaks  up  and 
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being  the  two  steering  lines,  is  virtually  intact.  As  a  result, 
fluids  particles  are  kept  outside  of  the  jet  by  the  outer  stable 
and  unstable  manifolds  which  form  barriers  to  mixing.  At  the 
same  time,  if  the  wave  speed  of  the  secondary  mode  is  smaller 
than  that  of  the  primary  mode,  this  boundary  breaks  up  and 
Lagrangian  particle  motions  can  penetrate  into  the  jet  via 
chaotic  mixing.  However,  even  in  this  second  case  the  inside 
of  the  jet  is  bounded  by  the  virtually  unbroken  inner  separa- 
trices  and  serves  as  a  barrier  to  transport  through  the  core  of 
the  jet. 

As  steering  lines  approach  each  other,  which  occurs  as 
one  goes  deeper  below  the  surface,  it  appears  that  mixing,  both 
within  and  across  the  jet,  should  be  facilitated.  Our  group  is 
studying  this  in  the  context  of  the  cusped  jets.  For  the  jet  with 
a  single  cusp,  the  steering  lines  approach  each  other,  but  the 
limit  is  quite  degenerate  as  they  will  converge  exactly  at  the 
cusp.  The  jet  clearly  becomes  very  thin,  but  indications  are  that 
it  remains  a  barrier,  as  seen  in  Figure  8. 

In  the  Bickley  jet  case,  ^  the  phenomenon  of  separatrix 
reconnection  occurs.  In  this  scenario,  heteroclinic  orbits  exist 
that  stretch  across  the  jet  between  the  two  rows  of  cat’s  eyes 
that  lie  either  side  of  the  jet  centerline.  This  would  appear  to 
be  a  level  at  which  mixing  is  strongly  facilitated  and  can  be 
realized  across  the  entire  jet.  However,  as  pointed  out  by  Pratt 
(private  communication)  this  may  be  hard  to  realize  in  an 
actual  situation  because  of  the  strong  potential  vorticity  gradi¬ 
ents.  Indeed,  this  is  corroborated  by  the  single  cusped  jet  case. 

The  phenomenon  of  convergent  steering  lines  can  also  be 
studied  in  the  context  of  the  double  cusped  jet  of  Pratt  et  al.^. 
Here  they  will  converge  as  the  speed  of  the  primary  wave  is 
reduced  to  match  the  speed  at  the  center  of  the  jet,  which  lies 
in  a  third  region  of  constant  potential  vorticity  between  the  two 
outer  regions.  In  the  case  of  a  smooth  jet  profile,  the  derivative 
U\y)  is  necessarily  zero  at  the  speed  at  which  the  two  steering 
lines  converge.  We  have  performed  a  local  analysis  to  study 
the  possible  dynamics  in  such  a  situation.  For  simplicity,  we 
assume  that  the  converging  steering  lines  occur  at  y=0  and  that 
both  (t)(y)  and  0'(y)  are  non-zero  at  y=0.  We  introduce  the 
following  change  of  variables: 

C  =  x:-ct,  Ti=a"’‘^y, 

Retaining  just  the  leading  order  terms  of  the  vector  field  results 
in  the  following  scaled  equations, 

^  =  (0)  /  -  a'^  O'  (0)  cos  C 

(9) 

^  O  (0)  sin  -  e  a  ^2  sm(k2^-a'^  an). 

At  e  =  0  there  are  homoclinic  loops  separated  by  a  distance 
k^  =  271,  and  heteroclinic  orbits  form  separatrices  on  only  one 


side  of  the  jet.  Figure  9  shows  the  resulting  flow  with  one 
neutral  mode  superimposed  and  the  parameters  set  to  u'\0)=2, 
O(0)=i,  o'(0)=i,  and  a=  0.01.  A  significant  issue  is  whether  the 
break-up  of  heteroclinics  on  the  outside  of  the  jet  still  occurs. 
By  using  the  Melnikov  function  from  (8),  we  have  been  able 
to  show  that  it  will  occur  in  a  uniform  manner  all  the  way  down 
to  the  convergent  steering  line  point.  This  indicates  that  the 
region  of  the  convergent  steering  line  presents  no  more  of  an 
effective  barrier  than  at  higher  depths. 

In  the  numerical  simulation,  when  a  second  neutral  mode 
is  added  to  the  streamfunction  i.e.,  e  >  0  is  set  in  (9),  both  the 
heteroclinic  and  homoclinic  tangles  are  easily  observed,  as 
seen  in  Figure  10.  The  characteristic  shape  of  the  heteroclinic 
tangle  is  easily  observable  in  the  figure  and  it  shows  that  the 
break-up  of  the  heteroclinic  orbit  results  in  significant  mixing. 
Further  numerical  experiments  also  indicate  that  particles  do 
indeed  traverse  the  entire  width  of  the  jet.  More  intricate 
calculations  and  additional  numerical  exp)eriments  are  needed 
to  ascertain  whether  it  is  a  region  of  more  efficacious  mixing. 
Of  particular  significance  are  simulations  for  somewhat  larger 
amplitudes  of  the  secondary  mode,  as  the  dynamics  in  that  case 
is  not  amenable  to  perturbation  methods  like  the  Melnikov 
method.  Related  results  for  two-dimensional  maps  suggest 
that  for  such  larger  perturbations  most  barriers  to  transport  are 
destroyed  and  new  mixing  mechanisms  are  created. 
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Introduction 

The  subject  of  double-diffusive  convection  covers  the 
study  of  fluids  in  which  there  are  gradients  of  two  or  more 
properties  with  different  molecular  diffusivities.  Many  inter¬ 
esting  phenomena  have  been  discovered  in  the  context  of  the 
ocean,  where  heat  and  salt  provide  the  competing  properties 
and  hence  this  subject  is  also  known  as  thermohaline  or 
thermosolutal  convection.  An  overview  of  the  area  is  given  in 
the  recent  articles  by  Schmitt  (1994,  1995).  The  following 
paragraph  takes  snippets  from  them  in  order  to  give  a  quick 
glimpse  of  the  history. 

Over  a  century  before  Melvin  Stem  discovered  salt  fin¬ 
gers,  W.  Stanley  Jevons  performed  the  first  salt  finger  experi¬ 
ment  in  an  attempt  to  model  cirrus  clouds.  Remarkably,  he 
seemed  to  realize  that  a  more  rapid  diffusion  of  heat  relative 
to  solute  played  a  role  in  the  experiments.  However,  he  over¬ 
simplified  the  physics  and  incorrectly  assumed  that  the  “inter¬ 
filtration  of  minute,  thread-like  streams”  was  a  general  result 
of  superposing  heavy  fluid  over  light  fluid.  Interestingly,  Lord 
Rayleigh  became  aware  of  these  experiments  more  than  two 
decades  later.  He  repeated  the  Jevons’  experiments  at  the 
Cavendish  Laboratory  in  Cambridge  in  April,  1880.  The  re¬ 
sults  led  him  to  initiate  the  study  of  buoyancy  effects  in  fluids 
by  formulating  several  stability  problems  for  a  stratified,  but 
non-diffusive  fluid.  His  neglect  of  diffusion  meant  that  he 


missed  an  opportunity  to  discover  double-diffusive  convec¬ 
tion.  The  modem  study  of  double-diffusive  convection  began 
with  Melvin  Stem’s  article  on  “The  Salt  Fountain  and  Ther¬ 
mohaline  Convection”  in  1960.  Stommel  et  al  (1956)  had 
suggested  that  a  flow  (the  salt  fountain)  would  be  driven  in  a 
thermally-conducting  pipe,  but  it  was  Stern  who  realized  that 
the  two  order  of  magnitude  difference  in  heat  and  salt  diffusivi¬ 
ties  allowed  the  ocean  to  form  its  own  pipes.  These  later  came 
to  be  known  as  “salt  fingers”.  Stem  also  identified  the  potential 
for  the  oscillatory  instability  when  cold,  fresh  water  overlies 
warm,  salty  water  in  the  1960  paper.  This  “diffusive-convec¬ 
tion”  process  was  demonstrated  later  by  Turner  and  Stommel 
(1964). 

For  the  purpose  of  this  article,  we  will  focus  on  the 
situation  with  warm  salty  fluid  below  cold  fluid.  The  fluid  lies 
between  two  horizontal  walls  placed  a  distance  L  apart,  and  is 
of  infinite  extent  in  the  horizontal  (x-y)  plane.  The  warm  fluid 
lends  to  rise  and  the  salty  fluid,  being  heavier  than  fresh  fluid, 
tends  to  fall,  so  that  if  we  begin  with  the  liquid  at  rest,  then 
increasing  the  temperature  difference  will  eventually  desta¬ 
bilize  the  rest  state  and  yield  either  an  oscillatory  motion  or  a 
steady  convective  motion.  In  particular,  we  will  address  the 
oscillatory  onset  of  instability  (the  “diffusive-convection” 
problem)  and  examine  a  family  of  patterns  that  may  arise.  We 
wish  to  find  new  solutions  that  bifurcate  from  the  rest  state, 
and  determine  their  stability  to  perturbations. 
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Methodology 

Stability  of  Fluid  Motions 

We  begin  by  identifying  the  velocity,  pressure,  tempera- 
ture  and  solute  concentration  that  satisfy  the  governing  equa¬ 
tions  and  which  represent  a  “base  flow”.  The  velocity  field  for 
the  base  flow  that  we  will  investigate  in  the  double  diffusion 
problem  is  zero  (the  “rest  state”),  the  temperature  and  solute 
concentrations  depend  linearly  on  the  vertical  variable  z,  and 
the  pressure  field  is  a  quadratic  polynomial  in  z.  If  this  base 
flow  is  disturbed  slightly,  then  the  disturbance  may  decay,  or 
may  grow,  depending  on  whether  the  flow  is  stable  or  unstable. 
If  the  flow  is  unstable,  then  the  solution  will  change  from  the 
base  flow  to  something  else,  but  to  what?  We  will  give  a  partial 
answer  by  considering  small  disturbances,  and  how  the  distur¬ 
bances  that  grow  the  fastest  interact  with  each  other  at  the 
conditions  close  to  the  onset  of  instability  (the  “critical  situ¬ 
ation”).  We  therefore  express  the  total  solution  as  a  superpo¬ 
sition  of  the  base  flow  and  a  perturbation,  and  write  down  the 
equations  that  govern  the  perturbation.  In  the  following,  a 
solution  will  usually  mean  a  solution  for  the  perturbation. 

spatial  Periodicity 

The  governing  equations  for  the  perturbation  to  the  base 
flow  show  that  if  one  stands  at  the  origin  and  looks  in  any 
direction,  the  problem  looks  the  same.  Many  bifurcation  prob¬ 
lems  in  fluid  mechanics  involve  one  or  more  spatially  un¬ 
bounded  directions  and  a  continuum  of  modes.  A  full 
description  of  the  set  of  bifurcating  solutions  in  such  a  context 
is  a  rather  formidable  problem  that  has  not  been  solved. 
Therefore,  one  typically  imposes  some  spatial  periodicity  on 
the  problem  and  then  confines  attention  to  solutions  satisfying 
this  given  periodicity.  For  instance,  one  may  look  for  solutions 
that  are  periodic  with  respect  to  the  vectors 

=  and  a:  2  =  M'- (0,1,0). 

These  vectors  are  said  to  span  a  hexagonal  lattice  of  period  W. 
In  the  X  and  y  directions,  the  solution  is  then  doubly  periodic; 
that  is,  at  any  time  t,  the  temperature  at  the  position  x  =  (x,y,z) 
is  equal  to  the  value  at  x+n^  Xj+n^  X2,  for  any  integer  np  n2. 
The  lattice  obtained  from  this  double  periodicity  is  invariant 
under  the  symmetries  of  the  hexagon;that  is,  rotation  by  mul¬ 
tiples  of  60  degrees,  reflection  across  the  vectors  a^  defined  by 

471  471  ,  1  V3 

and  reflection  across  axes  perpendicular  to  the  a-. 

The  same  periodicity  condition  holds  for  the  other  un¬ 
knowns:  the  solute  concentration,  velocity  field  and  pressure 
field.  This  simplifies  the  problem  immensely  and  allows  us  to 


use  a  Fourier  series  expansion  of  the  solution.  We  next  de¬ 
scribe  a  method  that  can  be  applied  to  find  ordinary  differential 
equations  for  the  dominant  Fourier  amplitudes,  and  thereby 
find  families  of  spatially-periodic  structures  on  a  planar  lattice. 
Rolls  and  hexagons  are  two  of  the  solutions  that  are  periodic 
on  a  hexagonal  lattice,  and  often  observed  in  convection 
problems.  Other  lattices  that  may  be  considered  are  the  square 
and  rhombus  (Golubitsky,  Swift  and  Knobloch,  1984;  Dionne, 
Golubitsky,  Silber  and  Stewart,  1994). 

Center  Manifold  Theoiy  and 
Dynamical  Systems 

The  analysis  of  spatially  periodic  bifurcating  solutions 
and  their  stability  proceeds  in  two  steps.  First,  the  governing 
equations  are  reduced  by  projections  to  a  finite  number  of 
nonlinear  ordinary  differential  equations  (looss  and  Joseph, 
1980,  Chow  and  Hale,  1982).  These  are  evolution  equations. 
The  center  manifold  theorem  is  invoked  for  this  purpose. 
Roughly  speaking,  the  theorem  says  that  in  the  neighborhood 
of  criticality,  the  dynamics  is  governed  by  the  interactions 
among  the  finitely  many  critical  modes.  This  places  the  prob¬ 
lem  in  the  framework  of  dynamical  systems.  At  this  stage, 
much  of  the  details  of  the  original  partial  differential  equations 
and  boundary  conditions  now  appear  in  the  coefficients  of  the 
evolution  equations.  Thus,  seemingly  different  physical  prob- 


Figuro  1 . 

The  neutral  stability  curve  for  the  oscillatory  onset  (line)  and 
steady  onset  (dashed  line)  of  instability  for  the  stress-free  problem 
is  shown.  The  Lewis  number  is  x  <  1. 
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Figure  2. 

We  display  the  four  standing  wave  solutions:  standing  rolls, 
standing  hexagons,  standing  regular  triangles,  and  standing 
patchwork  quilt  These  are  contour  plots  for  an  amplitude  func¬ 
tion,  e.g.,  the  vertical  component  of  the  velocity.  (Photographs  for 
these  figures  were  taken  by  Prof  David  Wagner,  Department  of 
Mathematics,  University  of  Houston,  Texas.) 


Standing  Regular  Triangles  Standing  Patchwork  Quilt 


lems  will  be  seen  to  reduce  to  the  same  types  of  evolution 
equations.  Secondly,  the  existence  and  stability  of  the  bifur¬ 
cating  solutions  are  determined  by  analyzing  the  system  of 
ordinary  differential  equations.  This  general  methodology 
handles  a  large  class  of  stability  problems,  where  the  onset  of 
instability  is  due  to  just  a  finite  number  of  modes  with  the  rest 
of  the  modes  being  stable,  as  is  typical  for  viscous  flows. 


A  thermal  Rayleigh  number  R^,  a  salinity  Rayleigh  number 
R5,  a  Prandtl  number  P,  and  Lewis  number  t  are  defined  by: 

= ga  (9* ,  -  0'  0)  LV(k^v),  =  g  ^  1  -  S*  0)  LV(KrV), 

P  -  V/Kj,  %  = 

where  g  denotes  gravitational  acceleration,  and  a  dimension¬ 
less  measure  of  gravity  is  G=gL^/k^^. 

A  basic  solution  is 

v  =  0.  0=1 -z,  5=l-z.  p=Po-Gz  +  (Rj-R^)P(z-^), 

where  Po  is  a  constant.  Since  we  are  concerned  with  the 
stability  of  this  rest  slate  and  with  solutions  bifurcating  from 
it,  we  set  up  the  problem  for  the  perturbations0,  p  and  v  = 
(u,  V,  w)  to  this  solution.  The  equations  satisfied  by  these 
variables  are: 

the  transport  equations  for  the  temperature  and  solute  concen¬ 
tration 


@  -w-Ag  =  -(v  V)0, 

(2) 

j  -  vv  -  tAJ  =  -  (v  ■  V)  S', 

(3) 

the  Navier-Stokes  equation  with  the  Oberbeck-Boussinesq 
approximation 

v-PAv4-Vp  +  (R^PS  -  R^jPQ)  -  (V  ■  V)  V,  (4) 
and  incompressibility 

div  V  =  0.  (5) 

Values  for  the  temperature  and  solute  concentration  are  pre¬ 
scribed  at  the  walls: 


Governing  Equations 

A  fluid  is  assumed  to  lie  between  two  walls  placed  at  z* 
=  0  and  z*  =  L  in  the  (x*,  y*,  z*)  plane.  Asterisks  denote 
dimensional  variables.  The  upper  plate  is  kept  at  a  constant 
temperature  0"‘o  with  solute  concentration  S^q,  and  the  lower 
boundary  is  kept  at  a  higher  temperature  and  a  higher 
solute  concentration  0*o  <  0*i ,  8%*  <  S*,  At  temperature 
0  Q,  the  fluid  has  thermal  coefficient  of  volume  expansion  a, 
solute  coefficient  of  volume  expansion  |5,  thermal  diffusivity 
solute  diffusivity  k^,  viscosity  |i,  density  p^,  and  kinematic 
viscosity  v=p/Po. 

Dimensionless  variables  (without  asterisks)  are  as  fol¬ 
lows:  (x,  y,  z)  =  (x*,  y*,  z*)  /  L,  t  =  Kjt*/L\  v  =  0  = 

(0*-0*o)/(0*i-0  V,S=  (S*-S  V/(S*rS*o),  P  =  P*lV(PoK^t)- 


0  =  0,  3^  =  0,  at  z  =  0  and  z=l.  (6) 

The  velocity  boundary  condition  at  the  walls  is  the  “no-slip” 
condition: 


=  0, 


at 


z  =  0  and  z  =  1. 


(7) 


Alternatively,  one  may  use  the  “stress-free”  boundary  condi¬ 
tion  (zero  normal  velocity  and  zero  shear  stress): 


fs  du  dw  ^  dv  dw  ^  ^  , 

w  =  0,^-f  — =  0,  — +  — =  0,  at  2  =  0,1. 
dz  dx  dzdy 


(8) 


For  the  stress-free  problem,  the  eigenfunctions  are  available 
in  closed  form  in  terms  of  trigonometric  functions  (Baines  and 
Gill,  1969,  Huppert  and  Moore,  1976,  Turner,  1979,  Nagata 
and  Thomas,  1984,  Knobloch,  1985,  Deane  et  al,  1987, 
Knobloch  et  al,  1986),  and  therefore  the  weakly  nonlinear 
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Figure  3. 

We  display  the  three  travelling  wave  solutions.  The  upper  pattern 
shows  travelling  rolls.  The  lower  patterns  are  the  travelling  patch- 
work  quilt  (1)  and  the  travelling  patchwork  quilt  (2).  The  arrows 
indicate  the  direction  of  travel,  (photographs  by  David  Wagner) 


k. 

Travelling  Rolls 


analysis  can  be  done  explicitly.  This  explicit  representation  is 
not  available  for  the  case  of  no  slip,  and  a  computational 
approach  is  then  used  for  the  analysis.  The  answer  to  the  question 
of  which  type  of  velocity  boundary  condition  is  pertinent,  lies  in 
which  phenomenon  one  is  trying  to  describe.  The  use  of  the 
stress-free  condition  is  relevant  for  modeling  situations  where  the 
fluid  layer  is  not  bounded  by  walls  but  by  free  surfaces;  for 
example,  double  diffusive  actions  may  be  taking  place  in  a 
sandwiched  layer  at  some  depth  say,  in  the  ocean. 

The  set  of  unknowns,  v,p,  J  and  0  and  0,  is  denoted  by  <I>. 
Equations  (2)  -  (8)  can  be  written  in  the  schematic  form 


L(5>  =  N2(.<J>,0), 


Here  the  operator  L  incorporates  the  linear  terms,  while  N2 
stands  for  the  quadratic  terms. 

Linear  Stability  Analysis 

For  the  case  of  infinitesimally  small  perturbations,  we 
linearize  (9).  The  rotational  symmetry  of  the  problem  about 
the  vertical  axis  shows  that  we  can  ignore  the  y  dependence 


and  pose  the  linearized  problem  in  the  x-z  plane.  The  method 
of  separation  of  variables  is  used  to  express  the  solutions  as 
exp(/ax+a()  multiplied  by  a  function  of  z.  The  real  linear 
operator  L  has  the  form  /i  B,  so  that  we  are  seeking  eigen¬ 
functions  C  which  satisfy  (A+  a  B)C  =  0.  The  wavenumber  of 
the  perturbation  is  a  and  the  growth  rate  is  Re  a.  We  wish  to 
map  out  the  regions  of  stability  (Re  rj  <  0)  and  instability,  by 
locating  the  neutral  stability  curve  (Re  a  =  0). 

Figure  1  shows  the  neutral  stability  curve  for  the  stress- 
free  problem  (Nield,  1967,  Baines  and  Gill,  1969)  for  the  case 
t  <  1,  which  is  typical  for  many  applications,  including  the 
ocean  media.  The  neutral  stability  curve  consists  of  two  lines 
meeting  at  a  Takens-Bogdanov  point.  The  steady  onset 
(dashed  line)  is  along 


277t^ 
4  ’ 


and  the  oscillatory  onset  is  along 


Rj-  — 


Tp  +1^ 


^Ps  +  (i+'')(l+p) 


277t‘* 

4 


The  critical  wavenumber  is  a—  %  /  T2,  the  same  along  the  entire 
neutral  stability  curve. 

For  the  no  slip  problem,  the  onset  conditions  need  to  be 
computed  numerically.  For  our  numerical  results,  we  use  the 
Chebyshev-tau  method,  which  approximates  the  eigenvalues 
o  with  infinite-order  accuracy.  This  is  a  spectral  method 
(Orszag,  1971,  Gottlieb  and  Orszag,  1983)  in  which  the  z-de- 
pendence  of  the  eigenfunctions  is  expressed  in  terms  of  Che- 
byshev  polynomials.  In  this  manner,  the  linear  stability 
problem  is  discretized  into  a  matrix  eigenvalue  problem.  The 
neutral  stability  curve  for  the  steady  onset  is 


^7.=^  + 1707.765,  a  =  3.12. 

The  neutral  stability  curve  for  the  oscillatory  onset  is  more 
complicated.  This  curve  lies  above  that  of  figure  1,  in  that  the 
critical  thermal  Rayleigh  number  is  higher  at  each  Rj.  More¬ 
over,  the  critical  wavenumber  increases  with  Rg. 


Nonlinear  Analysis 

Background 

In  the  two-dimensional  case,  there  is  a  critical  eigenfunc¬ 
tion  representing  waves  traveling  to  the  right  (with  wavenum¬ 
ber  a)  and  there  is  another  one  traveling  to  the  left  (with 
wavenumber  -a).  These  interact  nonlinearly  to  produce  trav¬ 
eling  and  standing  wave  solutions.  In  this  situation,  if  both  the 
standing  and  traveling  waves  bifurcate  supercritically,  then 
one  of  them  is  stable  (Ruelle,  1973).  Past  work  on  nonlinear 
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Figure  4. 

We  display  the  oscillating  triangles.  Upper  patterns:  t=0,  t-T/24, 
T  is  the  period.  Middle  patterns:  t=T/12.  t=T/8.  Lower  pattern: 
t^T/6.  (photographs  by  David  Wagner) 


t=0  t=3/24 


1=1/6 


analysis  has  been  on  the  stress-free  problem  (Knobloch,  1985). 
Deane  et  al  (1987)  derive  the  amplitude  equations  for  the  left- 
and  right-traveling  waves.  At  the  least,  they  need  to  calculate 
up  to  third  order  in  the  amplitude  functions.  The  calculation  of 
the  coefficients  in  the  equation  yields  the  Landau  coefficient 
which  determines  whether  the  solution  is  supercritical  or  not. 
They  find  that  the  third-order  Landau  coefficient  for  the  trav¬ 
eling  wave  turns  out  to  be  zero.  This  degeneracy  appears  to  be 
a  coincidental  one  (Knobloch  et  al,  1986)  due  to  the  adoption 
of  the  stress-free  condition.  The  standing  wave  solution  is 
supercritical,  indeterminate  and  subcritical,  respectively,  for 
Rg  less  than,  equal  to  and  greater  than  14585.4  for  the  case 
^-10-1/2.  p  =  At  the  fifth  order,  Deane  et  al  find  that  the 
traveling  wave  solution  bifurcates  supercritically,  so  that  there 
is  a  range  of  Rg  for  which  the  traveling  wave  solution  is  stable. 


What  are  the  third  order  Landau  coefficients,  and  what  is 
the  stable  solution  for  the  no  slip  case?  Computations  at  the 
above  values  of  Lewis  and  Prandtl  numbers  for  Rs=10^’ 
14585.4,  and  10^  reveal  that  the  solutions  are  both  unstable 
(Renardy,  1993).  As  expected,  the  selection  mechanism  is 
sensitive  to  the  boundary  conditions. 

What  other  periodic  solutions  are  there  in  the  fully  three- 
dimensional  planform  problem,  and  are  they  stable  to  the  more 
general  doubly  periodic  disturbances  in  the  x-y  plane?  The 
two-dimensional  solutions  we  have  discussed  correspond  to 
standing  rolls  and  traveling  rolls  in  three  dimensions.  The 
three-dimensional  stress-free  problem  has  been  studied  by 
Nagata  and  Thomas  (1984).  For  oscillatory  onset,  they  present 
the  linear  stability  of  standing  rolls,  standing  squares  and 
standing  hexagons  to  perturbations  restricted  to  have  the  same 
symmetry  as  the  solutions.  For  example,  they  examine 
whether  hexagon  pattern  solutions  are  locally  asymptotically 
stable  with  respect  to  hexagon  pattern  disturbances,  and 
whether  the  roll  solutions  are  stable  with  respect  to  roll-like 
disturbances.  The  next  issue  is  the  competition,  say  between 
rolls  and  hexagons,  and  the  investigation  of  the  traveling  wave 
solutions. 

Center  Manifold  Theorem  and 
Dynamical  Systems 

The  fluid  properties  are  regarded  as  fixed  and  we  define 
a  bifurcation  parameter  X=Ry-Rjc,  where  R-j^  is  the  thermal 
Rayleigh  number  and  R^^-  is  its  critical  value.  Solutions  O  with 
the  periodicity  of  the  hexagonal  lattice  are  sought.  Details  of 
this  calculation  are  found  in  Renardy  (1993). 

At  criticality  (^=0),  there  is  a  pair  of  complex  conjugate 
eigenvalues  ±  io)  of  the  linearized  problem,  and  each  of  them 
has  six  eigenfunctions.  For  X  near  0,  'ii(^)  denotes  the  eigen¬ 
value  which  arises  from  perturbing  -ico.  We  denote  by  Ck(^)» 
k=l,2,...,6  the  eigenfunctions  belonging  to  -li(X),  i.e. 

L(-4(?t))W  =  0, 

and  those  belonging  to  -  jl(X)  are  their  complex  conjugates.  We 
will  need  to  compute  the  eigenfunctions  at  X=0. 

The  critical  eigenfunction  computed  in  the  (x,z)-plane  has 
the  form  =  ^(z)  exp  (ia^  •  x),  where  aj=(a,0,0),  x  =(x,y,z), 
and  a  is  the  critical  wavenumber.  We  may  visualize  the  six 
eigenfunctions  as  follows.  The  vectors  ±  a^  ±  a2  and  ±  a3 
defined  by  equation  (1)  emanate  from  the  center  of  a  hexagon 
and  terminate  at  its  six  vertices.  The  critical  eigenfunctions  are 
waves  propagating  in  the  directions  of  a^,  a2,  a3  and  -  a^,  -  a2, 

-  a3.  We  denote  the  (complex)  amplitude  of  the  wave  propa¬ 
gating  in  the  direction  of  a^  by  z-  and  the  amplitude  of  the  wave 
propagating  in  the  direction  of  -  a-  by  z-^3. 

We  can  decompose  a  solution  O  in  the  form 

0  =  0i+% 
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(10) 


Figure  5. 


6 

Oi  =  2ReXz*Cyfc. 

1 

where  Zj  are  the  complex-valued  amplitude  functions  and  'P 
represents  a  superposition  of  eigenvectors  (and  possibly  gen¬ 
eralized  eigenvectors)  belonging  to  stable  eigenvalues.To  use 
the  Center  Manifold  Theorem,  we  need  to  have  the  fact  that  ^ 
consists  of  stable  modes.  The  expression  on  the  right  hand  side 
of  (10)  converges  to  O  in  the  mean  square.  The  ^  and  the 
functions  in  'F  form  a  complete  set.  We  wish  to  find  equations 
that  govern  the  evolution  of  the  amplitude  functions  z^.  In  order 
to  do  that,  we  need  to  find  an  explicit  approximation  to  'F  for 
small  perturbations;  at  this  stage,  we  invoke  the  Center  Mani¬ 
fold  Theorem. 

The  Center  Manifold  Theorem  states  that,  in  a  neighbor¬ 
hood  of  0=0,  there  is  a  manifold  T  (called  the  center  manifold) 
of  the  form  'F=t(Zi,Z2,Z3,Z4,Zg,Z6)  with  the  following  proper¬ 
ties: 

1.  All  solutions  with  initial  data  on  the  center  manifold 
remain  on  the  center  manifold  as  long  as  they  remain 
small. 

2.  All  small  periodic  solutions  lie  on  the  center  manifold. 

3.  The  stability  of  a  small  periodic  solution  is  determined 
by  its  stability  within  the  center  manifold,  in  other 
words,  all  Floquet  multipliers  corresponding  to  direc¬ 
tions  outside  the  center  manifold  are  stable. 

These  properties  imply  that  only  quadratic  terms  in  the  asymp¬ 
totic  approximation  to  the  center  manifold  are  required.  Thus, 
the  first  property  above  is  used  to  write  the  following  expres¬ 
sion 

6 

T  =  'I'2  + 'F2  =  2  Re  (  X  2/  z,-  Xy), 

ij=  1 

where  the  dots  indicate  terms  of  higher  than  quadratic  order. 
The  Yij  and  Xij  contribute  to  second  harmonics  and  the  mean 
flow  mode. 

The  second  property  above  says  that  the  solutions  of 
interest  (the  small  periodic  solutions,  taking  into  account  the 
nonlinear  terms)  are  found  on  the  center  manifold.  These  are 
solutions  close  to  the  linear  eigenfunctions.  There  need  not  be 
a  center  manifold  if  amplitudes  become  large. 

By  forming  appropriate  inner  products  and  projections 
onto  appropriate  spaces,  we  are  able  to  calculate  the  functions 
Vij  Xij‘  The  total  solution  O,  accurate  to  second  order,  is 
now  substituted  back  into  the  governing  equations,  leaving  the 
z-  as  the  only  unknowns.  After  further  projections,  using  the 
symmetries  of  the  hexagonal  lattice  and  the  theory  of  normal 
forms,  the  following  reduced  system  is  obtained: 


Wavy  rolls  (1).  Top  first  row;  t=0,  t=T/16,  T  is  a  period.  Second 
row  patterns;  t=T/8,  t=3T/16.  Third  row;  t=T/4.  t=5T/16.  Fourth 
row;  t=3T/8,  t=7T/16.  Fifth  row;  t=T/2.  (photographs  by  David 
Wagner) 


1=1/2 
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Figure  6. 

Twisted  patchwork  quilt  Upper  patterns:  t=0,  t=T/24,  Tisaperiod. 
Middle  patterns:  t=T/12,  t=T/8.  Lower  pattern:  t=T/6.  (photo¬ 
graphs  by  David  Wagner) 


t=l/6 


dzi 

—  +F,-  (Zi_Z2Z3_Z4,Z5,Z6,A,)  =0.J  =  (11) 


Fj  (zj^2'23,z4^5^g,X)  =  n-  (\)  z,  +  a,  (X)  I  z,  P  Zj  +  Oj  (X)  I Z4  P  Zj 

-KX3  (X)  ( I Z2 1^  + 1 Z3  P)  zj  +  a4  (X) 

(  I  Z5  P  +  I  Zg  P)  Zi  +  05  (X)  (Z2  Z5  +  Z3  Zg)  Z4  +  . . . 

The  dots  denote  terms  of  higher  order,  and  ai,...,a5  denote 
Landau  coefficients.  The  F-,  i=2,..,6  are  related  to  Fj  by  prop¬ 
erties  of  the  hexagonal  symmetry.  We  have  now  reduced  the 
governing  equations  to  a  system  of  six  complex  nonlinear 
ordinary  differential  equations. 


Pattern  Selection  on  the  Hexagonal 
Lattice 

In  (11),  we  set  u  =  and  obtain 

^  +  F(tt,X)  =  0,  (12) 

where  mgIR^^  and  X  is  a  real  parameter.  Starting  with  an 
equation  of  this  general  form,  Roberts,  Swift  and  Wagner 
(1986)  use  group  theory  to  classify  new  solutions  and  deter¬ 
mine  their  stability  properties  by  considering  their  symme¬ 
tries.  This  approach  has  allowed  them  to  find  a  larger  family 
of  solutions  than  was  previously  available.  Their  results  have 
a  generality  in  that  they  apply  to  other  problems  that  can  be 
expressed  in  the  form  of  (12)  and  have  the  same  symmetries. 
They  find  eleven  qualitatively  different  classes  of  bifurcated 
solutions:  (i)  Standing  rolls,  (ii)  Standing  hexagons,  (hi) 
Standing  regular  triangles,  (iv)  Standing  patchwork  quilt,  (v) 
Travelling  rolls,  (vi)  Travelling  patchwork  quilt  (1),  which  is 
always  unstable,  (vii)  Travelling  patchwork  quilt  (2),  (viii) 
Oscillating  triangles,  (ix)  Wavy  rolls  (1),  (x)  Twisted  patch- 
work  quilt,  (xi)  Wavy  rolls  (2).  These  are  solutions  which  have 
enough  symmetries  so  that  the  twelve  degrees  of  freedom  are 
reduced  to  two.  For  example,  invariance  under  rotations  by 
multiples  of  60  degrees  and  reflections  across  the  axes  of  the 
hexagon  yields  a  two-dimensional  solution  set  where  the  zi  are 
equal,  called  the  standing  hexagon.  The  conditions  for  stability 
of  each  solution  with  respect  to  perturbations  with  the  double 
periodicity  is  also  provided  in  Roberts  et  al.  To  apply  these 
results  to  a  specific  problem,  it  is  necessary  to  calculate  the 
coefficients  that  appear  in  (11). 

Figures  2  -  7  are  contour  plots  of  the  vertical  component 
of  the  velocity  for  each  of  the  solutions.  Solutions  (i)  -  (iv)  are 
standing  waves,  i.e.,  the  patterns  simply  oscillate  with  time, 
and  are  shown  in  figure  2.  These  plots  can  be  thought  of  as 
pictures  taken  at  time  zero.  Solutions  (v)  -  (vii)  are  traveling 
waves,  i.e.,  the  pattern  moves  to  the  right  without  changing  its 
shape,  and  are  shown  in  figure  3.  The  four  remaining  solutions, 
shown  in  figures4  -  7,  are  neither  standing  nor  traveling  waves. 
The  time  dependence  is  indicated  by  a  series  of  pictures.  For 
example,  figure  4  shows  the  oscillating  triangle  at  t=0,  1/24, 
1/12, 1/8  and  1/6  of  a  period.  The  pattern  evolves  through  the 
rest  of  the  period  before  going  back  to  the  t=0  picture. 

Numerical  Results  for  the  Oscillatory 
Onset  for  Double  Diffusion 

For  Lewis  number  T=10''^  and  Prandtl  number  P=  1,  we 
have  seen  that  traveling  waves  are  stable  in  the  two-dimen¬ 
sional  stress-free  problem  for  a  range  of  solutal  Rayleigh 
numbers  Rg,  while  it  is  unstable  for  the  no  slip  problem.  We 
pursue  the  no  slip  problem  as  a  model  for  the  fluid  bounded 
by  walls,  and  examine  the  stability  of  the  11  solutions  in  the 
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three-dimensional  case.  To  achieve  an  oscillatory  onset,  the 
solutal  Rayleigh  number  should  be  larger  than  approximately 
10^  results  in  the  range  10^<  Rs<10^  will  be  presented.  It  is 
found  that  the  solutions  are  unstable  until  the  solutal  Rayleigh 
number  becomes  relatively  large.  As  the  Rg  increases,  so  does 
the  Rj  at  criticality,  since  the  solute  concentration  effect  must 
counteract  the  thermal  effect  to  provide  for  the  oscillatory 
onset  (at  Rs=2  x  10^,  the  cridcal  R^  is  of  order  1  million).  For 
solutal  Rayleigh  numbers  between  8x10^  and  8.5  x  10^,  the 
standing  rolls  are  stable  and  the  other  10  solutions  are  unstable. 
There  is  a  window  of  Rg  from  9.0  x  10^  to  about  2.0  x  10^, 
where  all  11  solutions  are  unstable.  For  Rg«  4.0  xlO^,  either 
the  standing  hexagons  or  the  standing  regular  triangles  is 
stable  but  not  both,  and  the  remaining  9  solutions  are  unstable. 
As  a  reminder  of  the  scales  involved,  the  typical  properdes  for 
water  (Drazin  and  Reid,  1981)  are:  a  ~  5  x  lO'^^K'^  Oj  -  0^ 
<  10  K,  Kj  -  10'^  cmVsec,  v  -  10'^  cmVsec,  g«  10^  cm/sec^. 
Substituting  these  into  the  definition  for  the  thermal  Rayleigh 
number  to  be  1  million,  we  find  2,  L  being  measured  in 
cm. 

The  heating  of  cold  salty  water  is  more  appropriately 
modelled  by  a  Prandtl  number  of  P=10.0  and  Lewis  number 
x=0.01.  Numerical  results  were  obtained  for  solutal  Rayleigh 
numbers  up  to  order  a  million.  At  low  Rayleigh  numbers,  the 
solutions  were  found  to  be  unstable.  For  Rg  between  10^  to 
10*^,  the  standing  rolls  are  stable  and  the  other  solutions  are 
unstable.  For  Rg«  2x10“^,  three  of  the  standing  wave  solutions 
are  found  to  be  stable.  These  are  the  standing  rolls,  the  standing 
patchwork  quilt,  and  either  the  standing  hexagons  or  the 
standing  regular  triangles.  If  this  situation  were  realized  in 
nature,  we  may  see  a  coexistence  of  these  solutions,  a  patch  of 
one  lying  next  to  another. 

We  have  obtained  numerical  results  for  two  sets  of  Prandtl 
and  Lewis  numbers;  the  1 1  solutions  tend  to  be  unstable  for 
low  Rayleigh  numbers,  while  a  number  of  standing  wave 
solutions  may  be  stable  at  higher  Rayleigh  numbers. 

Outlook  for  the  Future 

The  center  manifold  reduction  method  is  applicable  to 
onsets  of  instability,  whether  steady  or  oscillatory  or  both,  in 
many  different  viscous  flow  situations,  and  continues  to  be 
applied  to  a  variety  of  problems.  The  group  theoretic  approach 
for  finding  classes  of  solutions  to  the  nonlinear  ordinary 
differential  equations  by  taking  advantage  of  the  symmetries 
in  the  problem  is  particularly  useful  for  convection  problems; 
recent  examples  include  ferrofluids  (Silber  and  Knobloch, 
1991),  viscoelastic  liquids  (Renardy  and  Renardy,  1992),  in¬ 
terfacial  stability  in  a  two-layer  system  (Renardy  and  Renardy, 
1988;  Joseph  and  Renardy,  1993)  and  magnetoconvection 
(Clune  and  Knobloch,  1994).  On  the  one  hand,  these  analytical 
results  serve  to  explain  observed  phenomena,  and  on  the  other, 
they  encourage  future  search  for  these  patterns  in  the  labora¬ 


tory.  For  instance,  we  have  found  that  three  different  patterns 
can  coexist  in  certain  parameter  ranges.  It  would  be  interesting 
to  see  which  one(s)  of  these  can  be  set  up  under  experimental 
conditions. 

This  article  has  taken  the  oscillatory  onset  for  the  double 
diffusion  problem  and  highlighted  the  dynamical  systems 
approach  coupled  with  group  theory.  Apart  from  the  oscilla¬ 
tory  onset,  the  double  diffusion  problem  can  have  a  steady 
onset  or  the  case  where  the  oscillatory  onset  merges  into  a 
steady  onset.  In  particular,  the  case  of  the  steady  onset  is 
relevant  for  the  investigation  of  salt  fingers.  The  subject  of  salt 
fingers,  and  their  stability,  is  an  ongoing  area  of  research, 
where  mathematical  modeling  together  with  the  type  of  analy¬ 
sis  presented  here  provide  fertile  grounds  for  future  interdis¬ 
ciplinary  collaboration.  For  instance,  why  are  square  patterns 
commonly  observed  (Proctor  and  Holyer,  1986)?  and  why  are 


Figure  7. 

lA'avy  rolls  (2).  Upper  patterns:  t=0,  t-T/24,  T  is  a  period. 
Middle  patterns:  t=T/12,  t=T/8.  Lower  pattern:  t=T/6.  (photo¬ 
graphs  by  David  Wagner) 


t=l/6 
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there  triangular  and  asymmetric  patterns  in  certain  regions  of 
the  ocean  (Schmitt,  1994b)? 
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